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ABSTRACT: 
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Provost 


The  two-  and  three-dimensional  unsteady  Navier-Stokes  equations  are 
solved  numerically  for  the  flow  field  about  an  impulsively  started  flat 
plate.   In  attempting  to  obtain  an  exact  time  dependent  solution,  several 
significant  results  were  observed.  First,  with  regard  to  the  formulation 
of  the  differential  equations  themselves,  it  appears  that  Poisson's 
equation  for  the  pressure  field  is  a  fundamental  equation  in  as  much  as  it 
allows  us  to  solve  for  pressure  most  exactly  at  any  given  time.  Secondly, 
the  difference  equations  must  be  carefully  and  consistently  formulated. 
In  this  research,  a  non -uniform  lateral  grid,  a  unique  interpretation  of 
the  continuity  equation,  and  "leap-frog"  integration  in  time  proved  to  be 
valuable  techniques  in  obtaining  an  exact  solution. 


LIST  OF  FIGURES 

Figure  No.  Page 

1.  Calculation  Region  for  Flow  Over  a  Finite  Flat  Plate 21 

2.  Time -Dependent  Solution  for  Flow  Over  a  Finite  Flat  Plate  ...   22 


li 


TABLE  OF  CONTENTS 

Section  Page 

NOMENCLATURE iv 

INTRODUCTION    1 

DIFFERENTIAL  EQUATION  FORMULATION      2 

DIFFERENTIAL  EQUATION  APPLICATION      8 

FINITE  DIFFERENCE  EQUATION  FORMULATION 10 

FINITE  DIFFERENCE  EQUATION  APPLICATION    11 

COMPUTER  PROGRAM Ik 

RESULTS  AND  DISCUSSION 16 

REFERENCES 19 

APPENDICES 23 

INITIAL  DISTRIBUTION  LIST Ik 


in 


NOMENCLATURE 

At  =  Finite  difference  time  step 

Ax  =  Finite  difference  grid  spacing  in  x 

Ay  =  Finite  difference  grid  spacing  in  y 

Az  =  Finite  difference  grid  spacing  in  z 

T)  =  Transformed  lateral  coordinate  =  e 

f  =  Arbitrary  dummy  dependent  variable 

L  =  Characteristic  length 

p,  =  Coefficient  of  molecular  viscosity 

v  =  Kinematic  viscosity 

p  =  Pressure 

p  =  Density 

Re  =  Reynolds  number  =  pVL/p, 

t  =  Time 

u  =  Velocity  in  x-direction 

v  =  Velocity  in  y-direction 

V  =  Characteristic  velocity 

w  =  Velocity  in  z-direction 

x  =  Axial  coordinate 

X  =  Body  force  vector 

y  =  Lateral  coordinate 

z  =  Cross -streamwise  coordinate 

u)  =  Successive  over -relaxation  parameter 

Subscripts 

i,j  =  Indicates  components  of  a  vector 
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NUMERICAL  EXPERIMENTS  WITH  THE  TWO -AND  THREE-DIMENSIONAL 

UNSTEADY  NAVIER  STOKES  EQUATIONS 

Gustave  Hokenson 

INTRODUCTION 

Since  Crocco  (Ref.  l)  suggested  seeking  asymptotic  solutions  to 
approximate  forms  of  the  parabolic  unsteady  Navier-Stokes  equations,  in 
lieu  of  attacking  the  elliptic  steady  equations,  many  researchers  (Refs.  2, 
3,  and  h)   have  successfully  applied  the  technique  to  solve  problems  ranging 
from  real  gas  nozzle  flow  to  the  hypersonic  blunt  body  problem.   The 
attempts  generally  are  aimed  at  retrieving  only  the  steady  solution  and 
this  limited  goal  allows  for  extensive  approximation  within  the  differ- 
ential and  difference  equations  as  long  as  the  steady  flow  is  exactly 
represented  by  the  equations.  This  approximation  of  the  unsteady  equations 
facilitates  the  numerical  solution  while  supposedly  not  interfering  with 
the  uniqueness  of  the  asymptotic  solution.  This  technique,  however, 
obviously  negates  the  intermediate  solutions  from  being  actual  representa- 
tions of  the  time  dependent  flow  field. 

The  essence  of  the  numerical  technique  is  to  specify  an  "arbitrary" 
initial  flow  field,  calculate  the  appropriate  spacial  derivatives  in  the 
Navier-Stokes  equations,  equate  their  sum  to  the  time  derivative  of  the 
respective  dependent  variable  and  step  forward  in  time.   Results  of  this 
investigation  revealed  that  the  combination  of  "arbitrary"  initial  condi- 
tions with  the  equation  approximation  previously  mentioned  creates  some 
numerical  problems  which  must  be  cloaked  in  various  numerical  subterfuges 
such  as  filters  and  arbitrary  intermediate  specifications  of  the  velocity 
field.   This  is  not  meant  to  slight  the  technique  for  it  does  provide  an 
easy  method  of  obtaining  apparently  unique  solutions  for  the  steady 


flowfield.   However  it  was  discovered  in  this  research  that  in  order  to 
guarantee  an  accurate  time  dependent  solution,  instantaneously  consistent 
with  the  equations  and  boundary  conditions,  much  more  care  is  needed  both 
in  the  formulation  and  the  numerical  implementation  of  the  differential 
and  difference  equations. 

The  purpose  of  this  investigation  was  to  obtain  the  exact  time 
dependent  solution  for  the  flow  about  a  finite,  infinitely  thin  flat  plate 
impulsively  started  in  its  own  plane  into  a  uniform  incompressible  flow. 
Both  the  two-  and  three-dimensional  Navier -Stokes  equations  were  applied 
to  the  problem  and  the  difficulties  characteristic  of  both  approaches  are 
presented.   A  variety  of  techniques,  all  consistent  with  the  equations  and 
the  boundary  conditions,  were  employed  to  obtain  exact  time  dependent  solu- 
tions without  recourse  to  non-rigorous  numerical  manipulations. 

DIFFERENTIAL  EQUATION  FORMULATION 
The  Navier -Stokes  equation  for  the  flow  of  an  incompressible, 
Newtonian  fluid  can  be  written  in  Cartesian  tensor  notation: 
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With  the  aid  of  Equation  1,  Equation  2  may  be  written  in  conservative  form 
which  is  most  amenable  to  exact  numerical  formulation  (Ref.  5).   In  conser- 
vative form,  Equation  2  becomes 
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Equations  k   and  1  form  a  closed  system  for  the  solution  of  the 
dependent  variables  u.  and  p  since,  with  the  appropriate  boundary  and 
initial  conditions,  there  are  as  many  equations  as  there  are  unknowns. 
These  equations  have  been  applied  directly  (Ref.  6)  to  obtain  steady  solu- 
tions to  the  Navier-Stokes  equations  with  one  of  the  momentum  equations 
serving  as  an  equation  for  pressure.  Early  in  the  course  of  the  investi- 
gation a  variety  of  problems  were  discovered  with  this  approach  simply  in 
retrieving  steady  solutions  and  its  utilization  to  obtain  the  desired  time 
dependent  solutions  proved  completely  unfeasible.   In  lieu  of  solving  a 
momentum  equation  for  pressure ,  Equation  h   may  be  differentiated  and  summed 
to  obtain  a  time  independent  Poisson's  equation  for  pressure  in  terms  of 
the  instantaneous  flow  field.   If  we  differentiate  Equation  k   with  respect 
to  x.  and  apply  the  Einstein  summing  convention  we  get 
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Assuming  the  functions  to  be  continuous,  we  can  interchange  the  order  of 
differentiation  in  the  temporal  convective  and  diffusive  terms  and  obtain 
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for  gravitational  body  forces  we  can  write 


a*i 


and  for   incompressible   flow,  Equation  1  states 
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Applying  these   conditions  to  Equation  6,  we  obtain  the  desired  Poisson's 
equation  for  pressure 
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With  Equations  1,  h,   and  9  we  have  a  system  of  equations  greater  in  number 
than  the  unknowns  and  we  must  choose  the  most  fundamental  set  to  solve. 
The  specification  of  the  apparent  optimum  choice  was  one  of  the  goals  of 
this  research  and  is  discussed  in  the  next  section. 

For  generality  in  problem  solving,  we  non-dimensionalize  Equations  1, 
U,  and  9  and.  for  simplicity  allow  the  new  variables  to  have  the  same  symbol 
as  the  terms  in  the  original  equations.   It  is  understood  from  this  point 
on  that  only  the  non-dimensional  quantities  will  be  discussed.   The  variable 
transformation  is  defined  by 
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With  this  transformation  the  governing  equations  can  be  written 
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where 


For  future  use  these  equations  are  expanded  for  the  three-dimensional  case 
in  rectangular  Cartesian  coordinates  as  follows 
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For  the  two-dimensional  case  we   can  simplify  these   equations  by  dropping  all 
terms    involving  w  and  z  and  their  derivatives.      In  this   case,  the   equations 
become 
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Just  as  the  conservative  form  of  the  convective  derivative  is  that 
which  is  most  appropriate  for  numerical  integration,  the  diffusive  .terms 
also  can  be  reformulated  to  most  exactly  satisfy  the  equations  when  finite 
differenced  (Ref.  5).   The  reformulation  of  the  diffusive  terms  is  accom- 
plished by  the  following  substitution  in  the  three-dimensional  equations 
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and  for  the  two-dimensional   system 
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where  the   appropriate   continuity  equation  has  "been  differentiated  to  obtain 
Equations   23-27.      After  performing  these   substitutions   the   3-d  equations 
become 
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For  the  two-dimensional  case  all  w  and  z  terms  and  their  derivatives  may  be 
dropped  and  we  retain  the  following  equations. 
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DIFFERENTIAL  EQUATION  APPLICATION 
As  was  mentioned  earlier,  there  is  a  choice  with  regard  to  the 
differential  equations  which  can  be  used  to  solve  the  problem.  As  a  general 
rule,  it  was  found  necessary  for  both  the  continuity  equation  and  Poisson's 
equation  for  pressure  to  be  satisfied  at  each  time  step,  including  t  =  0, 
to  guarantee  an  accurate  time  dependent  solution.  Based  on  this  principle 
each  streamwise  velocity  component  was  calculated  from  its  respective 
momentum  equation,  the  cross -streamwise  velocity  from  the  continuing  equa- 
tion and  pressure  from  Poisson's  equation.   If  there  is  no  distinct  cross - 
streamwise  direction,  for  example  in  a  wake  or  cavity  recirculation  region, 
a  momentum  equation  is  solved  for  each  velocity  component,  Poisson's  equa- 
tion is  solved  for  the  pressure,  and  continuity  is  checked  and  guaranteed  by 
the  use  of  correction  factors  (Ref.  7)  in  Poisson's  equation  in  the  terms 

^  _du.-,     ,   „  .  au. 
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which  cannot  be  arbitrarily  set  equal  to  zero  if  continuity  is  not  being 
used  to  solve  for  one  of  the  dependent  variables. 

For  the  particular  problem  we  have  attacked  the  differential  equation 
application  is  as  follows 
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u  from  u  momentum 
2-d    v  from  continuity  (37) 

P  from  Poisson's 

u  from  u  momentum 

w  from  w  momentum 
3-d  (38) 

v  from  continuity 

P  from  Poisson's 

In  order  to  obtain  exact  time  dependent  solutions,  it  was  found 
necessary  to  specify  physically  meaningful  and  compatible  initial  condi- 
tions throughout  the  flow  field,  not  "arbitrary  initial  conditions"  as  is 
so  often  quoted  in  the  literature.  The  compatible  and  physically  meaning- 
ful initial  conditions  for  this  investigation  were  specified  as  follows 

u  -  0  on  the  surface  of  the  plate 

u  =  1.0  everywhere  except  on  the  plate 
2-d  (39) 

v  =  calculated  from  continuity 

P  =  calculated  from  Poisson's 

u  =  0  on  the  plate 

u  =  1.0  everywhere  except  on  the  plate 
3-d    v  =  calculated  from  2-d  continuity  (ko) 

w  =  calculated  from  3-d  continuity 
P  =  calculated  from  Poisson's 

There  may  be  some  question  about  the  approach  for  the  3-d  problem,  however 
no  inconsistent  results  were  found  using  this  approach.   It  appears  that  the 
calculation  of  w  from  3-d  continuity  after  calculating  v  from  2-d  continuity 


would  result  in  a  null  matrix  for  w.   That  this  is  not  the  case  is  observ- 
able from  the  way  in  which  the  differential  continuity  equations  were 
implemented,  using  locally  spacial  averaged  values  of  u,  v,  and  w.  This 
procedure  is  outlined  in  the  next  section  and  the  appropriate  difference 
equations  are  listed  in  the  Appendices. 

FINITE  DIFFERENCE  EQUATION  FORMULATION 
Based  on  the  approach  outlined  in  the  previous  section,  the  spacial 
partial  derivatives  of  the  dependent  variables  in  each  of  the  differential 
equations  were  finite  differenced  using  three  point  central  difference 
formulas .   Special  care  was  taken  to  avoid  known  inaccurate  difference 
models  (e.g.  the  Du-Fort-Frankel  representation  of  diffusive  terms)  when  we 
are  seeking  exact  time  dependent  solutions.   Substitution  of  the  appropri- 
ate difference  expressions  into  the  momentum  equations  allows  for  the 
calculation  of  the  instantaneous  temporal  derivatives  of  the  dependent 
variables  u  and  w  based  on  knowledge  of  the  current  dependent  variables . 
The  differencing  procedure  for  the  momentum  equations  is  outlined  in 
Appendix  A,  followed  by  a  listing  of  the  respective  calculations  of  the 
temporal  derivatives  as  they  appear  in  the  computer  program.   In  both  the 
two  and  three-dimensional  situations  Poisson's  equation  was  used  for  the 
calculation  of  pressure.   Numerical  experiments  were  conducted  using  the 
y-momentum  equation  for  the  solution  of  P,  however,  since  we  dealt  with  the 
exact  equations  the  time  derivative  of  v  must  be  formed  and  added  to  the 
appropriate  spacial  derivatives  in  order  to  calculate  the  partial  deriva- 
tive of  P  with  respect  to  y.  With  the  pressure  on  the  freestream  boundary 
known,  it  was  in  theory  possible  to  integrate  the  equation  inward  from  the 
outer  boundary.   However,  numerical  instabilities  arose  from  the  inability 
to  calculate  physically  compatible  initial  conditions  since  at  this  time  v 
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is  known  only  at  one  time.   In  addition,  this  method  of  calculating  P  has 
inherent  accuracy  limitations  and  this  study  revealed  the  necessity  of 
calculating  P  to  a  very  high  degree  of  accuracy,  including  the  initial 
conditions.  This  discovery  served  as  the  impetus  for  using  Poisson's  equa- 
tion which  was  solved  using  classical  successive  over relaxation  techniques 
with  all  derivatives  in  the  equation  being  represented  by  three  point 
central  difference  models.  The  finite  differencing  of  Poisson's  equation 
is  presented  in  detail  in  Appendices  B  and  C  along  with  a  listing  of  the 
encoded  form  as  it  appears  in  the  full  computer  program. 

The  most  crucial  finite  differencing  formulation  which  was  encountered 
was  that  of  representing  the  continuity  equation  by  a  formula  which  is  both 
accurate  and  physically  meaningful.   The  method  which  was  settled  upon  is 
equivalent  to  a  third  order  Runge-Kutta  approach  and  the  details  of  its 
formulation  along  with  the  computer  listing  is  presented  in  Appendices  B 
and  C. 

In  each  of  the  equations  which  were  finite  differenced  the  cross  - 
streamwise  partial  derivatives  were  formulated  with  the  option  of  using  a 
non-uniform  grid  spacing.  This  allows  for  the  most  efficient  use  of 
computer  memory  and  concentrates  the  calculation  grid  in  the  regions  where 
the  dependent  variables  vary  most  rapidly.  The  central  difference  formulas 
for  both  the  first  and  second  order  partial  derivatives  with  a  non-uniform 
grid  are  developed  and  presented  in  Appendix  A. 

FINITE  DIFFERENCE  EQUATION  APPLICATION 
The  finite  difference  equations  discussed  in  the  previous  section  were 
applied  to  the  solution  of  the  incompressible  flow  about  an  impulsively 
started  flat  plate  in  the  manner  outlined  in  Equations  37,  38,  39?  and  kO . 
The  physically  compatible  initial  condition  field  was  found  and  the 
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appropriate  equations  were  used  to  continue  the  velocity  field  downstream 
in  time  with  Poisson's  equation  being  solved  for  the  new  pressure  field. 
When  explicitly  updating  the  velocity  field  in  time,  a  stable  solution  is 
obtainable  for  a  limited  range  of  time  step  increments.  Karplus  (Ref.  8) 
obtained  numerically  stable  solutions  with  a  similar  system  of  equations 
with  a  restriction  for  the  time  step  which  can  be  written  in  non-dimensional 
variables 

< 

At  ^ 5^7-n =i rr         (M) 


Ax   Ay   Az   Re  V  .  2    A  2    .  2 
^  Ax    Ay    Az 

Equation  Ul  proved  to  be  a  completely  satisfactory  stability  criterion  in 
this  work. 

In  an  effort  to  obtain  the  most  exact  solution  (to  second  order  in  At) 
the  mechanism  of  using  effective  central  differencing  in  time  was  introduced 
and  is  included  as  an  option  in  the  computer  program.   In  general,  an 
independent  variable  f  can  be  calculated  at  a  future  time  from  a  second 
order  Taylor  series  expansion  in  time  as  follows 


2  2 

f(x,y,z,t+At)  =  f(x,y,z,t)  +  ||  (x,y,z,t)  (At)  +  3_|  (XjVjZ  jt)^|l- 

at 

Since  f(x,y,z,t)  is  known  and  2_  can  be  calculated  from  the  differential 

ot  p 

3  f 
equations,  it  only  remains  to  calculate  — =■.   For  the  Navier -Stokes  equa- 

of  St 

tions  |_  is  an  explicit  function  of  the  dependent  variables  and  their 
ot  o 

spacial  derivatives  and  it  therefore  follows  that  — =■  may  be  calculated  by 

at* 

differentiating  the  governing  equations  with  respect  to  time   and  inter- 
changing spacial  and  temporal  derivatives.      However  the   implementation  of 
these  equations   is   lengthy  and  the  need  to  solve   a  second  Poisson's   equation 
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(for  — -)  increased  the  computing  time  intolerably.  Because  of  these 
aX, 

factors ,  it  was  decided  to  attempt  to  obtain  the  same  second  order  accuracy 
using  the  basic  difference  equations  appropriate  to  a  first  order  approach. 
This  was  accomplished  by  calculating  the  dependent  variables  at  time  t  +  At/2 
with  forward  finite  differencing  and  using  these  values  to  calculate  the 
spacial  derivatives  in  the  Navier-Stokes  equations  at  that  time.  The 
dependent  variables  are  then  stepped  forward  from  time  t  to  time  t  +  At 
based  on  the  temporal  derivatives  evaluated  at  t  +  At/2.  The  implementa- 
tion of  this  approach  is  most  clearly  understood  from  the  flow  chart  of  the 
computer  program  presented  in  Appendices  D,  E.  The  inherent  susceptibility  to 
numerical  instabilities  which  this  "leap-frog"  technique  exhibits  should  be 
damped  by  the  strongly  dissipative  terms  in  the  equations  (Refs.  9>  10,  an(i 

H). 

Finally,  with  respect  to  obtaining  the  solution  throughout  the  flow 
field,  some  technique  must  be  specified  to  meet  the  boundary  conditions. 
Apart  from  setting  u=v=w=0on  the  surface  of  the  plate,  some  procedure 
is  necessary  for  setting  u,  v,  w,  and  p  on  the  boundaries  of  the  calculation 
region.  Two  methods  are  available  and  have  been  proven  to  be  essentially 
equivalent.  First,  the  values  of  the  dependent  variables  which  have  been 
calculated  in  the  interior  of  the  region  may  be  extrapolated  to  the 
boundaries  using  the  arguments  of  symmetry  and  uniformity  across  the 
boundaries  to  generate  the  appropriate  extrapolating  function.  Second,  the 
equations  themselves  may  be  applied  to  the  boundaries  if  special  considera- 
tions are  taken  which  give  expression  to  those  values  needed  outside  of  the 
calculation  region,  based  again  on  symmetry  and  uniformity  arguments.  The 
latter  method  was  chosen  as  the  preferable  from  the  standpoint  of  consistency 
of  accuracy  and  error  generation  throughout  the  flow  field. 
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COMPUTER  PROGRAM 

The  computer  programs  which  encoded  the  solution  of  the  finite 
difference  equations  in  Appendices  B  and  C  are  listed  in  Appendices  D  and  E, 
along  with  the  flow  charts  appropriate  to  each.  With  regard  to  programming 
the  system  of  finite  difference  equations  and  several  mundane  considerations 
such  as  the  memory  capacity  of  the  computer  and  the  accuracy  required  of  the 
system  vs.  what  the  equations  are  demanding  of  it  are  deserving  of  some 
mention.   First,  with  regard  to  the  3-d  program  (with  the  inherent  3-d 
arrays  of  dependent  variables),  is  the  question  of  computer  core  storage 
capability.   The  minimum  number  of  3-d  arrays  which  must  be  dimensioned  is 
equal  to  the  number  of  dependent  variables.  Because  the  exact  difference 
equations  are  being  used,  the  calculations  of  new  variables  must  be  made  in 
conjunction  with  the  storage  of  old  ones  and  it  is  apparently  not  possible 
to  utilize  only  the  minimum  number  of  arrays .   However  in  this  program  the 
temporal  derivatives  are  written  on  discs  immediately  after  being  calculated, 
each  as  a  dummy  scalar.  When  the  calculation  of  the  dependent  variables  at 
a  future  time  (either  t  +  At  or  t  +  At/2)  is  being  initiated,  the  appropriate 
discs  are  rewound  and  the  scalar  strings  are  read  in  sequence  to  provide  the 
temporal  derivative  of  the  proper  dependent  variable  at  that  point  in  space. 

In  addition  to  the  artifice  of  using  virtual  memory  in  the  form  of 
simple  disc  writing,  an  additional  tack  was  tried.  This  was  to  map  the 
infinite  space  in  y  into  a  finite  space  with  the  transformation 

71  =  e"y 

if  this   is   done,   the   appropriate  partial  derivatives  become 
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i.-Tii 

by  '  ST1 


A  =   .2  j£_  +  .  A 

ay2  ^  ^2  ^ 


and  the  lateral  distances  can  be  stretched  such  that  there  are  fewer  points 
laterally  but  there  is  denser  data  in  the  region  of  steepest  gradients. 
This  transformation  also  allows  the  external  boundary  conditions  to  be  set 
somewhat  more  exactly  and,  in  conjunction  with  the  non -uniform  grid  size  as 
established  for  y,  should  give  the  optimum  use  of  total  lateral  space  per 
point  of  computation.  Because  the  results  were  adequate  at  the  Reynolds 
numbers  which  we  were  studying,  the  transformed  coordinate  was  considered  an 
unnecessary  complication  for  this  problem.   It  is  felt  however  that  many  of 
the  computational  difficulties  associated  with  the  exact  numerical  solution 
of  lower  Reynolds  number  problems  may  be  alleviated  by  the  use  of  this 
artifice. 

With  regard  to  the  questions  about  the  accuracy  of  the  solutions,  two 
specific  areas  of  investigation  were  dealt  with  in  this  research.  First  was 
the  utility  of  computing  in  double  precision  (real*8)  and  second  was  the  use 
of  higher  than  second  order  approximations  to  the  spacial  partial  derivatives 
The  computations  were  all  performed  with  a  variety  of  grid  sizes  and  accurate 
solutions  were  obtained  with  spacial  and  temporal  grid  sizes  large  enough  to 
allow  significant  differences  to  be  calculated  in  single  precision  (real*^). 
If  the  non-uniform  lateral  grid  is  set  up  with  extremely  small  spacing  near 
the  plate,  the  ensuing  small  differences  may  force  the  use  of  double 
precision  to  retain  the  most  accurate  solution.   In  general  it  is  felt  that 
this  is  the  only  situation  which  should  force  one  to  use  double  precision. 
Experimentation  with  the  use  of  higher  order  approximations  to  the 
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spacial  partial  derivatives  was  carried  out  with  the  following  two  representa- 
tions of  the  spacial  derivatives 

£    -  f.n*  +  8  f.i+i  -  8  f.i-i  +  h-2 
By  | .  12  toy 

&     -   fi+2  +  16  Vl  -  3°  f.i  +  16  f.1-l  -  f.1-2 

2  I  2 

ay  y  12  Ay 

The  results  using  these  models  were  identical  to  those  using  the  second  order 
approximations  with  adequate  grid  size  and  the  additional  complication  of 
increasing  the  thickness  of  the  boundary  on  which  the  boundary  conditions  are 
to  be  set  could  be  avoided  by  using  the  standard  formulas.  Therefore,  neither 
the  use  of  double  precision  nor  the  implementation  of  "more  accurate"  finite 
difference  models  was  found  to  be  of  general  use. 

RESULTS  AND  DISCUSSION 

The  exact  unsteady  two  and  three-dimensional  Navier-Stokes  equations 
have  been  applied  to  the  solution  of  the  incompressible  flow  about  an 
impulsively  started  flat  plate.   Some  typical  results  are  presented  in 
Figure  2  which  gives  the  transient  development  of  both  the  downstream  and  up- 
stream flow  on  the  center line  and  in  the  plane  of  the  plate.  For  our  purpose, 
it  is  not  particularly  startling  that  we  can  solve  this  specific  problem  but , 
more  importantly,  that  we  can  numerically  solve  for  the  time  dependent  flow 
field  using  the  two  and  three-dimensional  Navier-Stokes  equations.   In  so 
doing  a  set  of  non-rigorous  rules  has  been  obtained  which  can  be  used  as 
guidelines  in  the  application  of  the  equations  toward  solution  of  more 
complex  problems . 

The  Navier-Stokes  equations  written  in  strict  conservative  form  and 
physical  variables  can  be  integrated  to  obtain  accurate  time  dependent 

16 


solutions  providing  care  is  taken  in  the  formulation  of  the  problem.  At 
all  times  during  the  calculation,  the  continuity  equation  and  Poisson's 
equation  for  pressure  must  be  satisfied  to  ensure  a  valid  solution  while 
each  strearawise  velocity  component  must  satisfy  the  appropriate  momentum 
equation.   In  addition,  the  initial  conditions  must  satisfy  the  continuity 
equation  and  Poisson's  equation  for  pressure  in  order  to  allow  a  "natural" 
flow  development.  The  boundary  conditions  can  be  satisfied  either  through 
application  of  the  equations  to  the  boundary  with  the  use  of  symmetry  and 
uniformity  arguments  or  through  the  use  of  polynomial  extrapolation  from  the 
interior  calculation  region. 

In  essence,  Poisson's  equation  for  pressure  has  been  introduced  as  a 
fundamental  equation  in  the  numerical  solution  of  the  full  equations 
because  of  the  need  for  generating  a  compatible  initial  pressure  field  with 
a  time  independent  equation  and  also  because  of  the  need  for  calculating  the 
pressure  field  to  an  extremely  accurate  degree.  This  accurate  pressure  field 
calculation  and  the  specification  of  physically  compatible  initial  conditions 
eliminated  the  need  for  all  numerical  artifices  such  as  filters  to  cover  such 
purely  numerical  phenomena  as  velocity  overshoot. 

Second  in  importance  to  the  adequate  solution  of  Poisson's  equation  is 
the  satisfaction  of  continuity  equation  at  each  point  in  time  and  space. 
Apparently,  only  one  finite  difference  formulation  of  this  equation,  that 
presented  in  the  Appendices,  is  entirely  adequate  in  obtaining  exact  time 
dependent  solutions  of  the  problem. 

The  governing  equations  were  finite  differenced  using  a  central 
difference  approximation  for  all  spacial  derivatives.  The  DuFort-Frankel 
formulation  of  the  diffusive  terms  was  avoided  in  order  to  guarantee  the 
most  exact  representation  of  the  spacial  derivatives  at  each  instant  in  time. 
In  order  to  obtain  apparent  central  differencing  in  time  the  dependent 
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variables  at  time  t  +  At   were  calculated  using  the  spacial  derivatives 
calculated  at  t  +  At/2  where  the  t  +  A/2   values  were  obtained  by  pure 
forward  differencing.   To  obtain  stable  magnitudes  for  these  time  steps,  the 
approximate  equation  given  by  Karplus  was  used  as  a  guide  once  accurate 
spacial  grid  sizes  were  chosen. 

Beyond  the  formulation  of  the  numerical  problem,  the  operational 
problem  of  obtaining  sufficient  computer  memory  space  was  the  most  difficult 
to  overcome.   For  both  the  two  and  three-dimensional  problems,  the  use  of 
disc  storage  was  applied  to  help  alleviate  the  problem.  As  is  clear  from 
the  flow  diagram  of  the  program,  the  intermediate  calculations  were  stored 
on  discs  as  dummy  scalars  and  re-read  as  matrices  when  the  dependent  vari- 
ables were  being  updated.  This  technique  allows  the  use  of  a  minimum  number 
of  dimensional  arrays.   In  addition,  the  use  of  a  non -uniform  grid  and 
mapping  the  calculation  region  into  a  finite  space  were  applied  to  help 
alleviate  the  space  problem  and  optimize  the  number  of  calculation  points  at 
those  places  where  they  are  needed. 

Several  further  investigations  are  now  underway  to  determine  the  effect 
of  the  initial  conditions  and  the  equation  approximation  on  the  uniqueness 
of  the  asympototic  solution.   In  addition  the  two  and  three-dimensional 
instability  problem  is  being  studied  from  the  standpoint  of  comparing  the 
effect  of  two  and  three-dimensional  disturbances  and  finite  vs.  infinitesimal 
disturbances  on  the  classical  instability  maps.   Finally  the  full  compressible 
problem  is  being  formulated  in  a  manner  identical  to  that  presented  here. 
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FIG.  2 
TIME  DEPENDENT  SOLUTION  FOR  THE  FLOW  OVER 
AN  IMPULSIVELY  STARTED  FINITE  FLAT   PLATE. 
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APPENDIX  A 


Finite  Difference  Formulas  for  Non-Uniform  Spacing 


1    F 


j-l 


".    Vl 


-►  y 


We  let  F  =  Ay  +  By  +  C  and  expand  around  y.  to  find  A,  B,  and  C, 

J 


F    =  A  Dely  (j)  +  B  Dely  (j)  +  C 


F.  =  C 
J 


F    =  A  Dely  (j-l)  +  B  Dely  (j-l)  +  C 


Solving  for  A  and  B  we  get: 


A  = 


^F...  -  P.   F.    -  F.  -. 

J+l    J    J-l    J  I 


-1-  J^-L      J  +   J -J-      J   I 

Dely  (J)  +  Dely  (J-l)  L  Dely  (J)    Dely  (J-l) J 


B  = 


Dely  (J)  Dely  (J-l) 
Dely  (J)  +  Dely  (J-l) 


_F.  n  -  F.    F  .  _  -  F. 
J+l    J  _   J-l    J 

Dely2  (J)   Dely2  (J-l) 


Now  that  A,  B,  and  C  is  known,  we  can  differentiate  F  at  y  =  y .  and  it  is 

J 


clear  that 


^  =  B     &=   2A 
*y      'By2 
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APPENDIX  B 


Finite  Difference  Formulation  of  3-d  Navier-Stokes  Equations 


I.      u  Momentum  Equation 
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A=(U( 1  +  1, J,K)**  2-U< 1-1, J,K)**2) /<2.0E0*DELX) 

61=11  (I  t  J+1»K)*V(  I  ,  J+l  ,K)/  (DELY(  J  )**2  ) 

B2  =  U(  I  ,  J,KI*V(I  ,  J,K)M  1.0E0/DELYC  J-l)** 2-1  .  OEO/DELY  (  J)  M  2  ) 

B3  =  U(  I,  Jr-1,K)*V(  I  , J-1,K)/DELY(J-1)**2 

B=DELY( J)»DELY( J-l )*  <  B1+B2-B3  )/ ( DELY ( J ) +DELY ( J-l ) ) 

C=(U(  I,  J,K+1)*W(  I,  J,K+1)-U(  I  ,  J, K-1)*W(I  ,  J»K-1)  )/(2.0E0*DELZ) 

D=(PHI < 1+1, J,K)-PHI ( 1-1, J,K)  )/<  2.0E0*DELX) 

E2A  =  V(  1+1 ,J+1 ,K) /DELY( J)**  2 

E2B=V(I+1,J,K)*(  1.0E0/DELY(  J-1M«2-1.0E0/DELY(  J)**2) 

E2C=V(I+1 ,J-1,K) /DELY( J-l )**2 

E2  =  DELY(  J)*DELYU~m < E2A+ E2B-E2 C ) /( DEL Y ( J ) +DELY ( J-l ) ) 

E1A=V( 1-1, J+1,K)/DELY(J  )**2 

C1B  =  V(I-1,  JtK)» (1  .  OEO/DELY  (  J-l  )  **  2-1  .0E0/  DELY  ( J  ) **2  ) 

E1C=V( I-ltJ-ltK)/DELY( J-l)** 2 

E1=DELY(J)*DELV( J-l)* (E 1A  + E 1B-E 1C )/ ( DEL Y ( J )+DEL Y( J- 1 ) ) 

E=(E2-E1)/(2.0E0*DELX) 

F=(W( I+lfJ,K+l)-W( 1+1, J,K-1)-W(I-1,J,K+1)+W(I-1,J,K-1))/ 

(4.0E0*DELX*DELZ) 

G1  =  U(  I  ,  J+l.K)  /DELY(  J) 

G2=U(I, J,K)*< 1.0EO/DELY(J-1)+1.0EO/DELY( J) ) 

G3=U(I . J-l, K)/ DELY (J-l) 

G=2.0E0*(G1-G2  +  G3) /(DELY<  J)+DELY( J-l) ) 

H=(U( It J,K+1)-2.0E0*LHI,J,K)+U<  I , J ,K- 1 ) ) /DELZ**2 

UTEE=-( A+B+C+D+REI*(E+F-G-H)) 
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II.      W  Momentum  Equation 
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A=(U(I+l,J,K)*W(I+l,JtK)-U(I-l,J,K)*W(I-l,J,K))/( 2.0E0*DELX) 
S1=V( I » J+1,K)*W(I , J+ltK)/DELY( J)**2 

B2=V( I,J,K)*W(I,JfK)*( 1.0E0/DELY(J-1)**2-1.0E0/DELY( J)** 2) 
B3  =  V(I , J-1.K)*W( I, J-l ,K)/DELY( J-l  )**2 
B=DELY( J-1)*DELY( J)t(Bl+B2-BZ)/(DELY( J)+DELY(J-1) ) 
C=(W(I,J»K+1)**2-W(I,J,K-1)**2)/(2.0E0*DELZ) 
D=(PHI ( I, J,K+1)-PHI (I, J,K-1) )/( 2.0E0*DELZ) 

E  =  (U(H-1,  J,K+1)-U<  1+1 » J,K-1 )-U(  I-1,J,K  +  1 )+U(I-l,  J,K-1) )/ 
1(4.0E0*DELX*DELZ) 
F2A=V(I ,J+1,K+1 )/DELY(J)**2 

F2B  =  V<  I  ,J,K+1)M 1.0E0/DELY( J-l ) *  *  2-1 .OEO/DELY ( J )**2) 
F2C=V( I,J-1,K+1)/DELY( J-l)** 2 

F2=DELY(J)*CELY( J-l)»  (  F2A+F2B-F2  C  )/ (  DEL  Y  (  J  ) +DELY ( J- 1 ) ) 
F1A=V(  I  » J+1,K-1)  /DELY(  J)**  2 

F1B=V(  I  ,  J,K-l)t(  1.0/DELY( J-1)**2-1.0/DELY( J)**2) 
F1C=V(I ,J-ltK-l) /DELYt J-l)**2 

F1=0ELY(J)*DELY(  J- 1 ) * ( F 1A  +  F 1B-F 1C) /( DELY (  J)+DELY(  J-l )  ) 
F=(F2-F1)/(2.0E0*DELZ) 

G  =  ( W( 1  +  1, J,K)-W(I-1,J,K) ) /(2.0E0*DELX) 
H1=W( I, J+1,K)/DELY( J) 

H2=W( I, JtK)*(1.0EO/DELY(J-l)+1.0EO/DELY( J  )  ) 
H3=W( I , J-1,K)/DELY( J-l) 

H=2.0E0*(H1-H2+H3)/(DELY(  J ) +DEL Y( J-l ) ) 
WTEE=-( A+B+C+D+REIME  +  F-G-H)  ) 
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III.      Poisson's  Equation  for   Pressure 


BCD  E 


2  2  2  222222  -  ? 

°_Z  +  U.  +  LZ        p  u    +  a  v      a  wg  +  r_a_  pauv  :     _a_  rsw]     ami  "j 
Bx2     By2     ,Z2--L3x2      ay2      sz2     la*  La7  J     asLayJ     axazjj 


2     2     2 

a  p   a  p  a  p 

After  finite  differencing  — -  +  — -  +  — -  ,  we  find  that  the  Laplacian  can 

ax      ay      ay 

be  written 


2     2     2 
2-|  +  ^-|  +  2_|  =  Al  +  A2  +  A3  -  P(l,J,K)*H 

ax   ay   az 


combining  equations  and  we  can  solve  for  P 

P(I,J,K)  =  (A  +  B  +  C  +  D  +  2.0  (E  +  F  +  G))/H 

where  A  =  Al  +  A2  +  A3  and  H  are  defined  in  the  subsequent  listing.  Since 
equation  is  not  an  explicit  solution  for  the  pressure  field,  we  must  employ 
some  iterative  process  such  as  successive  over  relaxation  which  is  used  here, 
With  this  approach,  equation  is  written 

P(I,J,K)    =    (1   -  U))    P(I,J,K)    +fJu(A   +B   +  C+D   +  2.0(E+F+  G))/H 


A1=(PHI(I+1,J,K)+PHI( 1-1, J, K))/DELX**2 
A2  = (PHI (I, J,K  +  1)+PHI(  ItJ.K-lH /DELZ**2 
A3=2.0c0* (PHI  (I»  J+1,K)/DELY(  JI+PHKI,  J-1,K)/DELY(  J-l)  )/ 
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1(DELY(  J  )+DELY(  J-l)  ) 
A= A1+A2  +A3 

B  =  (U(I  +  lt  JfK)**2-2.0EO*U(I  ,  J  ,  K)  n**2+U  ( 1-1 1  J  ,  K)**2  ) /DELX**2 
Cl=V(Ii J+lfK)**2/DELYCJ) 

C2=(V( I,J,K)**2)*(1.0E0/DELY<J-1)+1.0/DELY( J) ) 
C3=V( I , J-l ,K)**2/DELY< J-l) 
C=2.0E0MC1-C2+C3)/(DELY(  J)+DELY<  J-l)  ) 

D=(W(I  , J,K+1)*^2-2,0E0*W( I, J,K)**2+W< I » J » K- 1 )**2 ) /DELZ**2 
E2A=U(  1  +  1,  J+l»K)*V(H-lf  J+1,K)/DELY<  J)M~2 

E2B=U(I+1, J,K)>V(I+1, J,K)*(  1.0F0/DELY(  J- 1  )--*2-l  .  OE  G7DE  L  Y(  J  )  **  2) 
E2C=U(I+1, J-1,K)*VU+1 ,J-1 , K)/DELY(J-1)**2 
F2=DELY(J)*DELY(  J-lH  (  E  2A*E  2B-E  2C  )  /(  DELY  (  J)  +DELY  (  J-l )  ) 
E1A=U(  I-1,J  +  1,K)*V(I-1,J-H,K)/DELY(J)4?  2 

Eie=U(  I-lt  J,K)»V(I-lt  JtK)*  (l.OEO/ DELY  (J-l  )**-2-l  .0  EO/ DELY  ( J  )**2) 
E1C  =  U(  I-1,J-1,K)*V(  I-1,J-1,K)/DELY(  J-l)  **2 
E1=DELY(J)"*DELY< J-l)* ( E 1A+E 1B-E 1 C  )  / ( DEL Y( J  )+DELY( J-l)  ) 
E  =  (E2-E1)/(2.0E0*DELX) 

F2A=V(  I  ,  J  +  l,K  +  l)*W(If  J  +  1,K+1)/DELY(  JH  t  2 

F2B  =  V(I  ,J,K+1)*W<  I,  J,K  +  1)*(1.0E0/DELY(  J-1)*'*2-1.0E0/DELY(J)**2) 
F2C=V(I ,J-1,K+1)*W( I,J-1,K+1)/DELY( J-l)**2 
F2=DELY(J)*DELY( J-l  )M  F2A+F2B-F2C  )/ (  DEL  Y  ( J)+OELY( J-l)) 
F1A  =  V(  I  ,  J+lfK-1)  *W(I  f  J+1,K-1)/DELY(J)**2 

F1B=V(  I,J,K-1)*W(I,J,K-1)M 1.0E0/DELY(J-1)**2-1.0E0/DELY( J)* *2) 
F1C=V(  I  , J-ltK-l)*W(I, J-1,K-1)/DELY(J-1)**2 
F1  =  DELYU)*DELY(  J-l)* (F 1A+F 1B-F 1C) / ( DEL Y ( J ) +DELY ( J-l ) ) 
F=(F2-F1)/(2.0E0*DELZ) 

G2=U(I*-1,  J,K+L)*W(  I+lf  JfK  +  1  )+U(  I-l,Jf  <-l  )*W(I-1,  J,K-1) 
G1=U( I- 1, J,K+1)*W(I-1,J,K+1)+U( I+1,J,K-1)* W(I+l,JfK-l) 
G=(G2-G1)/(4.0E0*DELX*DELZ) 

H=(2.0EG7DELX**2+2.0E0/DELX**2+2.0E0/lDELY(J)*DELY(J-l) ) ) 
APHI=( 1.0E0-0MEGA)VPHI( I , J , K) * OMEGA* ( A+ B+C+D+2. OEO* (E+F+G) )/H 
RESID=ABS( APHI-PHI (I, J, K) ) 
ERR=MAX1(EPR ,RESID) 
PHI(I,J ,K)=APHI 
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IV.      3-d  Continuity  Equation  to  Solve   for  w 


Bw  _        pdu  +  Sv 
Bz  |_dx       By 


UEX2=U( I+l,J+ltK)+U(I+l, J-l, K)+U( 1+1, J+l, K+l)+U( 1+1, J-l, K+l) 
UEX1=U( I-1,J+1,K)+U( I-1,J-1.K)+U( I- 1T J+1,K+1)+U(I-1,J-1,K+1) 
UEX=( UEX2-UEX1)/(8.0E0*DELX) 

Vl=(V(I+l,J+l,K)+V(I+l,J+ltK+l)+V(I-lfJ+l,KJ+V(!-lfJ+l,K+l))/ 

1 4 . 0E0 

V2  =  ( [  V(  1+1, J,K)tV(Itl, J,K+1)+V(I-1,J,K)+V(I-1, JtK  +  1 J )/4.0E0 
V3=(V(  I+1,J-1,K)+V(  1+  1,  J-  l.K+U+VU-1,  J-1,K)  +  V(I-1,J-1  ,K+1)  )/ 
1 4 .0E0 
VY1=V1/DELY( J)**2 

VY2=V2* (1.0E0/DELY( J) **2-l. OEO/DELY( J-l)** 2) 
VY3=V3/PEIY(J-1)**2 

VtfY=DELY<  J)*DELY(  J-l)*(  VYH- VY2-VY3  )  /  (  DELY  <  J  J+DELY(  J-l)  ) 
W(  I,  J,K  +  1)=W(I,  J,K)-DELZ*(UEX+VWY) 


To  evaluate  r —  at   I,J,K  +  l/2  the  following  method  was  employed.      First  v 
dz  ' 

at  y.,-|  j  y.j  y.   -,   was   found  by  averaging  the  four  v  values   in  the  x-z  plane 

at   a  given  y.     With  these  average  values,   the  non-uniform  grid  derivative 

3v  / 

formula  was   applied  to  form  —  at   I,J,K  +  1/2.      In  a  similar  manner,   an 

average  u  at   I  +  1  and  I   -  1  was   found  by  averaging  the  four  values    in  the 

y-z  plane  at   each  I.      These  average  u's  were   then  central  differenced  to 


form  —  at   I,J,K  +  l/2.     Finally,   using  these  derivatives   evaluated  at 
dx 


Bu 

ax 

I,J,K  +  l/2,  the  w  field  is  calculated  using  Equation    using  a  central 
difference  method. 
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V.   2-d  Continuity  Used  to  Solve  for  v 


dv     3u 
By    "ax 


The  identical  procedure  was  applied  to  the  finite  difference  solution  of 
this  equation  in  two  dimensions  as  was  outlined  for  the  3-d  continuity 
equation.  The  principle  difference  being  that  the  appropriate  averaging 
of  the  dependent  variables  now  occurs  on  lines  instead  of  planes  in  the 
3-d  matrix  of  grid  points. 


C1  =  U(  I  +  1,J+1,K)+U(I  +  1,J,K) 
C2=U(I-1, J+l,K)+U(I-lt JtK) 
V(If  J+ltK)=V(Ii  J,K)-DELY(J)*(C1-C2)/K.0E0*QELX» 


VI.      3-d  Continuity  Used  to  Solve  for  v.     Following  the  procedure  out- 
lined in  Section  IV  of  this  Appendix,  we  can  interpret  the  three-dimensional 
continuity  equation  to  be  an  equation  for  v,  here  again  using  averaged 
values  of  the  dependent  variables   in  the  difference  formulation. 


UEX2=(U(I  +  lf  J,K-l)+U(H-lf  J+ltK-1)    +U(I  +  1»J»K+1)+U(I  +  1»J+1,K+1))/ 
14.0E0 

UEX1=(U(I-1,J,K-1)+U( I-lt J+ltK-l)+U(I-lf J,K+l)+U(I-lf J+l , K  +  l  )  )/ 
14.0E0 

UEX=( UEX2-UEX1)/(2.0E0*DELX) 

WZ2=(W(  H-1,J,K+1)+W(  I  +  lt  J+ltK+l)  +  WU-l,  J,K+1)+W(I-1,J+1  ,K+1)  )/ 
14.0E0 

WZ1  =  (  W(  I  +  l,Jt  K-l)+W(I+l,J+lfK-l)+W(I-lf  Jf  K-D+WC  I-1,J  +  1,K-1)  )/ 
14.0E0 

WZE=(WZ2-WZ1 )/(2.0E0*DELZ) 

V(  I,J+1,K)=V(I,  JTK)-DELY(J)<f  (UEX  +  WZE) 
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APPENDIX  C 
Finite  Differencing  of  2-d  Navier -Stokes  Equations 
I.      u  Momentum  Equation 

ABC  D  E 

2  ^2 


B_u 
St 


f  B_u_  +  3uv  +  B_P       2.  XjL  rSv~l       5  uV, 

L  dx         By        ax       Re    l3x  LByJ   "  .    2j  J 

oy 


A=(U( 1+1, J)**2-U( I -It  J )**2)/(2.0E0*DELX) 

BL=U(I t J+1)*V(I i J+l )/DELY( J )**2 

B2=U( I, J)*V( It J)*( 1.0E0/DELY( J- 1 )** 2-1. OEO/DELY ( J)**2) 

B3=U( I, J-l)*V(It J-1)/DEI_Y(J-1)**2 

B=DELY( J)*DELY( J-l) M B1+B2-B3) / ( DELY( J ) +DEL Y ( J-l ) ) 

C=(PHI(H-l,J)-PHI(I-ltJ))/( 2.0E0*DELX) 

D2A=V(I+lt J+1)/DELY(J)**2 

D2B=VMI+ltJ)*(  1.0E0/DELY(J-l)*i<2-1.0E0/DELY<J)**2) 

D2C=V(I+lt J-1)/DELY(J-1)**2 

D2=DELY( J)*0ELY( J-l)* ( D2A+D2B-D2C )/ ( DELY ( J ) +DELY ( J-l ) ) 

D1A=V(I-1,J+1)/DELY(J)**2 

D1B  =  V(I-1,  J)  +  (l  ,0E0/DELY(J-1)**2-1.0E0/DELY(J  )**2) 

D1C=V( 1-1, J-l) /DELY (J-l) **2 

D1=DELY(J)*DELY(  J-l)t  ( D1A+D 1B-D 1C ) /( DEL Y(  J)+DELY(  J-l)) 

D=(D2-D1)/(2.0E0*DELX) 

E1  =  U(  I»  J+1)/DELY(  J) 

E2=U(I, J)*( 1.0E0/DELY(J-1)+1.0E0/DELY( J) ) 

E3=U(I , J-1)/DELY(J-1) 

E=2.0E0*(E1-E2*E3)/(DELY( J)+DELY( J-l) ) 

UTEE=-( A+B+C+REI*(D-E ) ) 


II.      Poisson's  Equation  for  Pressure 


^2„       ^2^ 

.^2   2          ?2 

a  p  ,  a  p  , 

2            2 

[  B   u          B   v 
L       2               2 

Bx         By 

ax           By 

Bx  L  By  JJ 


2     2 
B  P  L  3  P 


From  finite  differencing  — =■  +  — p  ,  we  can  show  that 

Bx    By 
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2     2 
L_|  +  L_|  =  A1  +  A2  -  P(I,J)*H 
dx    dy 


Substituting  this  into  equation,  we  obtain 

P(I,J)  =  (A  +  B  +  C  +  2.0  D)/H 
To  obtain  a  form  appropriate  to  successive  over -relaxation,  we  can  write 
P(I,J)  =  (1  -  oo)  P(I,J)  +  »(A  +  B  +  C  +  2.0  *  D)/H 


A1=(PHI  (  1  +  1,  J)+PHI(  I-l.J)  )/DELX**2 

A2=2.0E0*(PHI  (I,  J  +  1)/DELY  (JI+PHK  I  f  J-l  I  /DEL  Y(  J-l )  )/ 
1(DELY( J)+DELY( J-l) ) 
A= A1+A2 

B=(U( 1+1, J)**2-2.0E0*U< I, J)**2+U(  I- It  J )**2) /DEL X**2 
C1=V( I , J+1)**2/DELY(J) 

C2=V( I, J)**2*( 1.0E0/DELY(J-1)*1.0E0/DELY( J) ) 
C3  =  V(I  , J-1)**2/DELY(J-1) 
C=2.0E0*(C1-C2+C3) /(DELY( J ) +  DE L Y( J-l ) ) 
D2A=U(I+1, J+l  )*V(  I+lt  J+1)/DELY(  J)**  2 

D2B=U(I+1, J)*V(I+1 , J)*(1.0E0/DELY<J-l)**2-i.0E0/DELY(J )**2) 
D2C=U( 1+1, J-1)*V( 1  +  1, J-1)/DELY( J-l)** 2 

D2=DELY(J)*DELY(  J-l)*  ( D2A+D2B-D2C )/ ( DEL Y ( J  )+DELY(  J-l)) 
D1A=U(  1-1  t  J+1)*V(I-1,  J-H)/DELY(  J)**2 

D1B=J(I-1,J)*V( 1-1, J)* ( 1.0E0/DELY(J-1)**2-1.0E0/DELY(J)**2) 
D1C=U(I-1, J-1)*V( I -It J-1)/DELY( J-l)**2 

D1=DELY(J)*DELY{ J-l ) M D1A+D 1B-D 1C ) /( DEL Y ( J) +DELY ( J-l ) ) 
D=(D2-D1)/(2.0E0*DELX) 

H=2.0E0/DELX**2+2.0E0/(DELY(J)*DELY( J-l ) ) 
APHI=(1.0E0-3MEGA)-*PHI(  I  ,  J)  +  0ME  GA*  ( A*B+C+2.  0E0* D  )  /H 
RESID=ABS(APHI-PHI (I, J ) ) 
EPR  =  VIAX1(ERR,RESID) 
PHKI,  J  )  =  APHI 


III.      2-d  Continuity,  Solving  for  v 

dv  _       du 
by  "    "  dx 

To  form  —  at  I, J  +  l/2,  average  values  of  u  at  I  +  1  and  I  -  1  were  found 
ay 

by  averaging  the  values  of  u  along  appropriate  J,  J  +  1  lines.  Finite 

Bu 
dx 


differencing  of  these  average  values  yields  —  at  I,   J  +  l/2. 


UEX2=(U(I+1»J)+U( 1+1, J+l) J/2.0E0 
UEX1=(U(I-1,J)+U( 1-1, J+l) )/2.0E0 
UEX=(UEX2-UEX1)/(2.0E0*DELX) 
V(I,J+1 )=V(It J)-DELY(J)#UEX 
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APPENDIX       D 
I.  FLOW  CHART  FOR  3-D   PROBLEM 


SET    I.   C. 

ITER=1 

0 

TEST  ICALC 

1 

f 

ITER=0 

STORE  U,  W 

4 

^ 

COMPUTE   Uf,  Wf 

ITER=1 

STORE  Ut,Wf 

M 

TEST  ITER 

0                 * 

w 

1 

DELT=DELT  1 

DELT=DELT  2 

0 

TESTICALC 

1 

r             4 

f 

1 

READ  U,   W 

READ    Uf,Wt 

0 

COMPUTE    U,  W 

COMPUTE    V 

COMPUTE    P 

TEST  ITER 

1 

^ 

r\\  iTniiT 

UU  I 

ru  i 

V 
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APPEND  I  X    D 

II.     LISTING    0-    COMPUTER    PROGRAM     FOR    3-0    CASE 

THE    FOLLOWING    FORTRAN    COOED    FINITE    DIFFERENCE    COMPUTER    PROGRAM 
SOLVES    THE    pULL    TIME    DEPENDENT     NAVIER    STOKES    EQUATIONS     FOR    THE 
THREE     DIMENSIONAL    FLOW    ABOUT    A     FINITE,     INFINITELY    THIN    FLAT    PLATE 
IMPULSIVELY    STARTED    IN    ITS     OwN    c>LANE.    "LEAP-FROG    TIME-WISE    INTEG- 
RATION     IS    INCLUDED    AS    AN    OPTION.     EXTENSIVE    DISC    WRITING    IS     EMPLOY- 
ED  TO    SAVE    COMPUTER    CORE    SPACE.     ADDITIONALLY,    THE     EQUATIONS 
THEMSELVES    ARE    APPLIED    TO    THF    BOUNDARIES     OF    THE    CALCULATION    REGION 
IN    LIEU    OF    EXTRAPOLATING    THE    VARIABLES    CALCULATED    WITHIN    THE 
REGION    TO    THE    BOUNDARIES. 

THE    FOLLOWING    PARAMETERS    MUST    Be     SPECIFIED 
ICALC    LE    ME    NF    LP1    LP2    NP  i    MP  2    RE    DELX    DELY    DELZ    OMEGA     ERRTOL    f>E     OF 


fc 

t 

rv 


ICALC  SPECIFES  THE  USE  OF  "LEAP-FROG"  IN  T  *■ 
LE  IS  THE  LENGTH  OF  THE  REGION  IN  DELX  STEPS  * 
ME     IS    the    HEIGHT    OF    THE    REGION     IN    DELY    STEPS    * 

*  # 

NE     IS    THE     WIDTH    OF    THE    REGION    IN    DELZ    STEPS 
LP1     IS    THE    START     OF    THE    PLATE 

LP2  IS  THE  END  OF  THE  PLATE 
NP1  IS  ONE  SIDE  OF  THE  PLATE 
NP2     IS    dNE    SIDE    OF    THE    PLATE 

ft 

RE  IS  THE  REYNOLDS  NUMBER  BASED  ON  U  AMD  L  \ 

DELX  IS  THE  GRID  SPACING  IN  X 

»• 

DELY    IS    THE    GRID    SPACING    IN    Y 

v. 

DELZ    IS    THE    GRID    SPACING    IN    Z 
OMEGA     IS    THE    RELAXATION    PARAMETER     IN    POI SSONS    ' 
ERRTCL    IS    THF    ERRCR    TOLERANCE     IN    PCI SSONS 
PE    IS    THE    FREESTREAM    PRESSURE     IN    PSF 
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c 

o 

f 

c 

c 

c 
c 


c 
c 
c 


QE  IS  THE  FREESTREAM  DYNAMIC  PRFSSORE  IN  PSF 

4        NOTE  THAT  NE  MUST  BE  AN  ODD  INTEGER 
'<■  r 

NCTE  THAT  OELY  IS  AN  ARRAY  DELY(ME) 

r-  -3 

ALL  VARIABLES  ARE  DIMENSIONED  F(LE,ME,NF)    * 

f  * 

>      ■<  •  +  ••-...;,*•■  *f  S*  ?    '  •  <"  <   *  *r  >  i»  $  >  >*  £  £  **»*.  £* .'«  -;    *■   • :  r  ••  s»«*  ..^.•^••.^•:i- 

DIMENSION    U(t0t30,2i), V(40t  30T  21),W(40t  30.21),DHI ( 40t3O?21) 
DIMENSION    DELYC3  0) 

Lf;i=LE-l 
LE2  =  LF-2 
LE3=LE-3 

MEJ.=I*E-1 
MF2=ME-2 
ME3=ME-3 

NE1  =  NE-1 
NE2-NE-2 

ME3=NE-3 

NPM=(NE+1) /2 

NPM1=NPM-1 

N»M2=NPM+1 

REI=1.0tO/RE 

PHIE  =  PE  /(2.0E0--QE) 

A=1.0EC/DELX 
B=1.0EO/DfcLY( 1) 
C=1.0EO/DELZ 

delt=i.qeo/(a  +  0.1eo-  b+g.1e0-c+2  .oeo*rei*  (  a**2+b**2+c*-<  2)  ) 
celti=delt 

DzLT2=DELTl/2.0EG 
IOO 

INITIALIZE    THE    CALCULATION    FIELD 
D  G    2     I  =  1  ,  L  E 
DC     2    J-1,ME 

DC    2    K=1,NE 


3U 


U(  If J,K  )  =  1.0E0 
V(I,J,K)=O.CEO 
W(I,J,K)=O.OEO 
PHKIf  J»K)=PHIE 
CONTINUE 

SFT    PLATE    TO    ZERO 
DC    3    I  =  LP1*LP2 
DC    3    K  =  NP1,NP2 
U(I,l,K)=O.OEO 
CCNTINUE 


V  EL  DC  IT  Y 


CALCULA 
DO  4  J  = 
DG  4 
DO  4 
C1  =  U( 
C2=U( I- 
V(I , J+l 
CGNTINU 


V    WITHIN    THE    REGION    FROM    2-D    CONTINUITY    EQ 


1  = 
K  = 
1  + 


TE     INITIAL 

lfMEl 

2, LEI 

I,  NE 

If J+1,K)+U(I+1, J,K) 

1, J+1,K)+U( 1-1, JfK) 

,K)=V(I , J,K)-DELY(J  )  -<C1-C2)/(4.GE0*DELX) 

F 


SET    UPSTREAM    V    HELD 

DO    5   J=1,ME 

DO    5    K=1,NE 

V(  l,J,K)=O.OEO 

CCNTINUE 

EXTRAPOLATE    DOWNSTREAM    V    FIELD 

DC    6    J=1,ME 

DC    6    K=1,NE 

V(LE,J,K)  =  3.0EO*  V(  LEI,  J,K)-3.0E0i'V(LE2,  J,K)+V(LE3,  J,K) 

CCNTINUE 

CALCULATE    W    FIELD    WITHIN    CALCULATION    REGION    FROM    3-D    CONTINUITY    EQ 

CC    7    1=2, LEI 

DO    7    J=2fMFl 

DO    7    K=NPM,NE1 

UEX2=U( 1+1 ,J+1,K)+U( I+1,J-1 ,K)+U( I+l» J+l» K+l)+U( 1+1, J- 1,K+1) 

UEX1=U( 1-1, J+lf K)+U(I-1,J-1 fK)+U< I-lf J+l f K+1)+U<I-1 » J-lf K+l) 

UEX=(UEX2-UEX1 )/ ( 8  .OEO-DEL X ) 

Vl=(V(I+lf J+1,K)+V(I+1 , J+l, K+l )+V(I-l,J+l,K)+V(I-l, J+l, K+l) )/ 
14. CEO 

V2=(V(I+1, J,K)+V( 1+1, J,K+1)+V( I-1,J,K)+V( 1-1, J, K+l) )/4.0FO 

V3=< V(  1  +  1  ,J-1 ,K)  +  V(I  +  1  ,J-1,K+1)+V(I-1,J-1 , K)+V(I-1, J-1,K+1) )/ 
14.0EO 

VY1=V1/DELY(JJ**2 

VY2=V2'( 1.0EO/DELY( J)**2-l. OEO/DELY( J-l  )**2  ) 

VY3=V3/DELY(J-1K*2 
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VwY=CFLY( J)  DELYt J-l) * ( VY 1+ VY2- VY3) /( CE L Y( J)+DELY( J-l)  ) 
w( I ,J,K+1) =W(I t J, K) -DELL*   (UEX+VwY) 
r  CCNTINUE 

CALCULATE  M  IN  THE  PLANE  OF  THE  PLATE 

DC    8    1  =  2,  LEI 

DC    3    K=NPM,NF1 

J  =  l 

UEX2=U(I  +  1,J  +  1,K)  +  U(I  +  1,J+1,K)  +  U(I  +  1,J<-1,K+1)+U(I+1,J+1  ,  K+l  ) 

UEX1=U(  1-1  ,J  +  1,K)+U(  I-1,J  +  1,K)+U(  1-1, J  +  1,K  +  1)+U(  I~l, J+1,K+1) 

UCX  =  (  UEX2-UEX1J  /(8.0E0*  DELX) 

Vl=(V(I  +  l,J  +  l,K)+V(I  +  l,J  +  l,K+l)  +  V(I-l,J+i,K)  +  VU-l,J+l,K+l))/ 
14. GEO 

V2=(V(I+l,J,K)+V(I+l,J,K+l)+V(I-l,J,K)+V(I-l,J,K+l))/4.0EO 

V3=(V(I+lt  J  +  lfK  J+V(I  +  l,J  +  l,K+l)+V(I-l,J+l,K)  +  V(I-lf  J+ltK+U)/ 
14. DEC 

V2=-V3 

VY1=V1/DELY(J )'  ^2 

VY2  =  V2-  ( l.OFO/DELYl J) *-  -2-1. OtJ/ DELYC J )     '  2) 

VY3=V3/DELY( J)C \  2 

VfoY=OcLY(  J)    OFLY  (  J  )  .  ( V Y 1+VY 2-V Y 3 ) / ( DEL Y ( J  >+DELY(  J  )  ) 

w( I,J,K+1) =W(I , J,K)-rELZ^ (UtX+VWY) 
:    CONTINUE 

CALCULATE    W    AT    THE    TGP    OF    THE    CALCULATION    REGION 

DC    9    1=2, LEI 

DC    9    K=NPM,NE1 

J=ME 

UtX2=U(  I+1,J,K)+U(  I  +  lt  J-1,KH-U(  I+1,J,K+1)+U(I+1,J-1,K+1) 

UrXl=U(  1-1  ,J,K)+U(  1-1  , J-l ,KT+U( I -1, J, K+l )+U(I-l, J-l, K  +  l ) 

UEX^(UEX2-UPXl)/(  8.0E0-  D^LX) 

Vi  =  (V(I+l,J,K)+\M!+l,J,K+l)+V(I-l,J,K)+V(I-l,J,K  +  l))/ 
14.0E-0 

V2=(V(I+1,J,K)+V(I+1,J,K+1)+V(I-1,J,K)+V(I-1,J,K+1))/4.0E0 

V3  =  (V(I+1,J-1,K)+V(I  +  1,J-1,K  +  1)+V(I-1,J-1,K)  +  VU-1,J-1,K+1))/ 
14.0EC 

VY1  =  V1/  DELY(J  )•'  -2 

VY2=V2<   (l.CEO/CELYU)  *    2-1  .OED/  DELY(  J-l  ;*•   2) 

VY3=V3/PELY(  J-l)'  -:  2 

VWY  =  LELY(  J)- *  DELY  (J-l)-*  (VYl  +  VY2-\/Y3)/(DFLY(J  )  +  DELY(  J-l)  ) 

W(l,  JfK  +  l)=W(I  »  J,K)-DELZMUF.X+VWY) 
•    CONTINUE 

INSERTION    OF     LECT    HALF    REGION    W    VALUES     BY    SYMMETRY 

DO    10     I  =2, LEI 

DC    10    J=2,ME1 

DC     10    K=NPM2,NE 

K1=NE+1-K 
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c 
c 

r 


rt(I,J,Kl)=W(I,J,K) 

10    CCNTINUE 

CALCULATE    W    CN    THE    FRCNT    OF 
DC    11    J=l, MF 


11 


13 

810 

150 

6  II 

14 


THE    CALCULATION    REGION 


CO     11    K=1,NC 
W(lt JtK)=O.OEO 
CCNTINUE 

CALCULATE    W    AT    THE     EXIT    3C    THE    CALCULATION    REGIJN 

DO    12    J=1,ME 

DC    12    K=1,NE 

w(LE,J*K)*3.0ED.^(LEl,  J,K)-3.0E0!' W(LE2t  J,K  )+W(LE3t  J»K) 

CONTINUE 


ON    THE     PLATE    SURFACE    EQUAL    TO    ZFRO 
I=LP1,LP2 


SET    W 

DO    13 

CO    13    K  =  NP1,NP2 

W(I ,  1,K) =0.0E0 

CONTINUE 


CALCULATE     INITIAL    P    FIELD 

WRITE(6,810) 

FOPMAT( 1H0,20X,' INTO    IC    PCI SSONS • ) 

GO    TO     200 

IC  =  1 

WRITE( 6,811) 

FORMAT ( 1H0,20X,  'OUT    CF     IC    POI SSONS') 

ITER=1 

CONTINUE 

IF(  ICALC.EQ.O)    GO    TO    18 

ITER=0 

REWIND    16 

REWIND    18 


C 

c 

c 
c 

c 

STORE  PREVIOUS  FIELD  ON  TAPES 
WRITE (16)  U 
WRITE  (18)  W 

18 

CONTINUE 

REWIND  19 

REWIND  20 

CALCULATION  WITHIN  THE  REGION 

CO  19  1=2, LEI 

DO  19  J=2,ME1 

DO  19  K  =  2fN€l 
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0'  (G1-G2  + 
, J, K+l )-2 
( A+B+C+D+ 
19,905)  J 
( 1PE14.7) 


I+lt 
1+1. 
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Y  (J) 


UT  F 
2-U( 
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t  JiK 
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)  /DE 
(  1.0 
)  /DE 
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1/PE 
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—  a  (  I 
Z) 

EL  Y( 
3  EC/ 
:  LY( 
33)/ 
.GEO 
RE  I 
TEE 


KCM    30    MOMENTUM    =  QUATinv 

I  - 1 ,  J ,  lO  >     2)/(2. GEO -DELX) 

+  1  f K)/(DELY( J  J  "     2) 

)-  (  l.OEO/CELYt J-l)        2-1.0E0/DELY( J) -     2) 

-ltK)/DELY  (J-l)       2 

(B1+B2-B3)/ (DELY(J)+DELY(J-i ) ) 
f<+l)-U(I,J,K-l).  W(  I  ,  JtK.-l)  )  /(  2.0L0    DELZ) 
(  1-1,  J,K  )  )/(2  .GEO4  DELX  ) 
L  Y  (  J  )  -;     2 

EO/DELY(  J-l)-1  2-1.0F0/DELY(  J)'      2) 
LY(  J-l  )"      2 

)     (E2A+E  2S-E2C) /  (  DE L Y ( J ) +DE LY ( J-l ) ) 
LY( J)>     2 

EG/DELY(  J-l)-*  ■  2-1.  0  EG/ C  ELY  (J  )* ■   -2  ) 
L  Y(  J-l)        2 
)     (  ElA  +  El&-elC)/"(  OEL  Y(  J  )+DELY(  J-l)  ) 

LX) 

+lf Jf K-l J-W( I-lfJ,K+l)+U(I-ltJtK-l)) / 

J) 

DELY(J-lH-1.0tO/DELY(  J  )  ) 

J-l  ) 

(DELY(  J)+DELY(  J-l)  ) 

J(  I, J,K)+U(  I, JtK-1)  */PELZ:      2 
(E+F-G-H) ) 


CULA 
U(  1  + 
V(  I  , 
V(  I, 
V(  I  , 
ELY( 
W(I  , 
PHI  ( 
U(  1  + 
GEO  ■ 
=  V(  I 

=  y/(l 
=v(  I 

GELY 
=  V(I 
=V(  I 
=V(  I 

DELY 
F2-F 
W(  I  + 


TION 
1  t  Jf 
J  +  l» 
Jf  K) 
J-l  t 
J-l) 
J,K  + 
If  Jt 
It  J, 
DELX 
t  J  +  i 
tJtK 
7  J-i 
(J) 
t  J+l 
tJtK 
t  J-l 
(J) 
l)/( 
It  J, 


OF 

*)  ;  W 

K)*  'f 
W  (I 

K)  »  W 
DEL 
1)^ 
K+l) 
K+l) 
DEL 
t  K  +  l 
+  1  ) 
t  K+l 
GELY 
tK-1 
-1) 
tK-1 
-ELY 
2.0~ 
K  )  -W 


WT 
(  1  + 
(  It 
t  Jt 
(I  t 
Y(  J 
2-W 
-PH 
-J( 
Z) 
)  /P 
(1  . 

)  /c 
( J- 

)/D 
(  1. 
)/D 
(  J- 

0    D 

(  I- 


FROM 
1  t  J  t 
Jf  1, 
K)  »  ( 
J-lt 
)  (E 
(  I,  J 
!  (J  t 
I+lt 

FLY( 
OcO/ 
ELY( 
1)  ( 
ELY  ( 
C/DE 
ELY( 

i  )*  ( 

ELZ  ) 
It  Jt 


3D    MOMENTUM    EQUATION 
K)-U(  1-1  t  J,K)    W(!-l,  JtK)  )/  (2.  DEC-   DELX) 
K)/D5LY( J) >  <  2 

1.  CEO/ DELY  (J-lp  >  2-1.  CEO/DEL  Y(  J)  >   '2) 
K)/DELY( J-l) > -2 
1+B2-B3) /( DELY( J)+DELY( J-l) ) 
,K-1  )■•*    2)/  (2.0E0-DELZ  ) 
J,K-1) )/( 2.0EG-DELZ) 
J,K-l)-U(I-l,JtK+l)+U(I-lfJtK-l))/ 

J)-     2 

PELY( J-l)'     2-1.0E0/DELY(J) •      2) 

J-l) •  «2 

F2A  +  F  28-F2C) /(DELY( J)+DELY( J-l) ) 

JM   >2 

LY(  J-l)«  '  2-1.0/DELY(  J)  f'2  ) 

J-l)  --2 

F1A+F18-F1C)/(DELY(J)+DELY(J-1)) 

K  )  )/(  2. GEO'  DELX) 
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c 


c 
c 


CALCUL 
DO    20 
DC    20 
K  =  l 

CALCUL 
K-l  TE 
A  =  (U(  I 
B1  =  U(  I 
32=U(I 
B3  =  U(  I 
B=CELY 
C  =  (U(I 
D=(PHI 
E2A=V( 
E2B=V( 
E2C=V( 
E2=DEL 
E1A=V( 
E16=V( 
E1C=V( 
E1=DEL 
E=(E2- 
FM  W(  I 
1(4. OEO 
G1=U( I 
G2  =  U(  I 
G3-U(  I 
G=2*0E 
H=(U( I 
UTEE=- 
WRITE  ( 
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-B3)/ 
(  I,  J, 
K)  )/( 
2 

(J-i) 
—  2 
E2B-E 
2 
(J-l) 


TUM    EOUATICN 

UNIFORMITY    ARGUMENT 
/  (2. OEO .-  D<=|_X) 
J)**  2) 

Y(J-1  )    ■  2-1.0E0/DELY( J)'     2) 
-i)~- 2 

( DELY( J)+DELY( J-l) ) 
KKW(I»  J,K)  )/(2.0E0*DELZ) 
2.0E0*DELX) 


>  "■  2-1.0E0/CELY(J-)*"2  ) 
2C)/(DFL  Y(  J  )+DELY(  J-D) 
-    2-1.0E0/DELY(  J)*>  21 
1C)/(DELY( J)+DELY( J-l) ) 
-W(I-1 , J,K+1 )+W (1-1, J, K) )/ 


ElB-E 


LY(J-1)  +  1.0E0/DELY( J)  ) 
1) 

FLY  (J  )+DEL 
( r,J,K)+U( 
+  F-G-H)  ) 


Y(J-l)  ) 

I  ,  J,K))  /DELZ<     2 


CALCULATION    OF    WT     FROM.    3D    MOMENTUM    EQUATION 

A  =  (U(  1  +  1,  J,K)  *w(  I  +  1,J,K)-U(I-1  ,  J,K)'  W(I-1  ,  J,K)  )/  (2  .OEO-  DELX) 

Bl*V(If J+l,K)*W(If J+ltK)/DELY(J )**2 

B2  =  V(I  , J,K) *W(I, J,K)^ '  (l.OEO/ DELY  (J-l)*  '•<  2-  1  .OEO/ DEL Y( J ) *  -2) 

B2=V(  I  ,  J-l,K)t  W(  I  ,  J-1,K)  /DELYl  J-l)«.  :  2 

B=CELY( J-l )    DELY (J )" (Bl+B2-B3)/( DELY( J)+DELY(J-1) ) 

C  =  (W(I  ,  J,K  +  1)«  <  2-W(  I,  J,K)  *«  2)/<2.0E0~DELZ) 
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D  =  (PH 
E=(U( 
U.OE 
F2A  =  V 

FZE=V 
F2C=\, 
F2  =  Dr 
F1A  =  V 
Flb  =  V 
F1C  =  V 
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p  =  (^2 
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OEO/QELY( J) + l.OEO/DELYC J) J 

ELY(J) 

33) /(DELY( J)+DELY(J) J 

.  OrO    U(  It  J,K)+tl{  I  tJtK)  )  /DELZ       2 

REI ■  (E+F-G-H)  ) 
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B2- 

r-     5   = 
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FIB 
F  K 
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hl  = 
H2= 
H3  = 
H  =  2 
WTF 
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CULA 
U(I  + 
V(I  » 
V(  If 
VU  , 
-H3 
ELY( 
W(  I  , 
PHI  ( 
U(  1  + 

oeo- 

=  V(I 
=  V(I 
=  V<  I 
=  -F2 
DELY 
=  V(  I 
=  V(I 
=  V(  I 
=  -Fl 
DELY 
F2-^ 
W(  1  + 
W(I  t 
W<  I» 
w(  I  1 

.OEC 
c  — ( 
TE(2 

TIMU 


TICN    CF    WT    FROM    3D     -iC'itNTUM    EQUATION 
ltJtK)*rf{  1  +  1,  J,K)-U(  i-L  ,  J.K)  'W(  I-l,J,lO)/(2.0i 
J+1,K) - W(I  t J+l, K) /DBLY( J)     -2 


2  L  X ) 


J,K)'  *U,J,!<)     (  1.0=0/D"rLY(  J) 
J+1,K)     M(  I  ,  J+l  t  !<)/.)  ELY  (J  ) 


2-1.0E0/Di.LY<  J)       2) 


J)OELY(J )     (61+B2-B3)/(l6LY(J)+DBLY( J) ) 

J,K+1M     2-w(I t J, K)     ^2)/(2.0FC^0ELZ) 

It J,K+1)-PHI< I, J,K) ) /( 2. OEO    DELZ) 

itJ,K+l)-U(I+lfJfK)-U(I-lfJfK+l)+U(I-l,JfK))/ 

DELX    DELZ) 

,  J  +  1»\+1)/DELY(  J  )  •  ■  2 

T  Jf  K+l)  -■  (1.0EO/DELY(  J)     -  2-1 .0 FG/ DEL Y < J  )     -  2) 

,  J+l, K+l)  /DELY(  J)  '  *2 

C 

(J)     CELY(  J)+(F2A+F2B-F2C)/  ( DELY ( J ) +CELY I J ) ) 

, J+1,K)/2ELY( J)        2 

,J,K)  -  (1  ,0/DELY  (J)--     2- 1.0/  DEL  Y(  J  )     *  2  ) 

,J+lfK)/D2LY( JH A2 

r 

Tj)  -  PELY  (  J)-   CF1a+FIB-F1CI/C0cLY{ J J+PELYiJ ) ) 

l)/(  2.03.0-  DELZ) 

l,J,K)-MI-l,J,K))/(  2.  GEO- DELX) 

J+l  ,KJ  /DELY (J) 

J,K J     (  l.OEO/DELYl J  )  +  1.  OF  O/DL  l_Y(  J)  J 

J+l.KJ/DELY (J) 

> (H1-H2+H3) /(DELY( J) +DE LY(JJ) 

A+B+C+D+REI •  (E  +  F^G-H)  ) 

0,904)     *TEE 


CALCULA 

Dr.  25  I 
J=l 

K^NE 
J-l    AMD 
C^LCuLA 
A=(  U(  1  + 

ei=u(  i » 

32-U( I  f 
B3=U(  I  , 
b3=-B3 

b=DELV{ 
C=(u<  I, 
D=(PmI  ( 
E2A  =  v'(  1 


=2,LC1 


A    CORNER 


K  +  l 
"ICN 


TE: 


MS    HAVE    BEEN    ADJUSTED 

UT     FRCM    3D    MnwiENTUM    EQUATION 

1,  J,K)**-2-U(I-l,  J,K) -    2) /(2.0E0    DELX) 

J  +  1,K)'  VUtJ  +  ltKj/(OFLY(  J)--   2) 

J,K)     V(I  »  J,K)     (1  .OEO/DLLY  (J)  -    2-1  .0£0/D: 

J+ltK)'  V(I  ,  J+ltM  /OcLYC  J)     *2 


LY(J>       2) 


J)  *DELY(  J)"    (B1+d2-B3  )7  (  DELY  ( J  ) +DELY  ( J  )  ) 

J» K)     W( I , Jf K)-U(I t J,K-1)      WU  .J.K-1)  )  /(2.0E0* DELZ) 

I+1,J,K)-PHI(T-1,J,K))/(2.0E0*<.2LX) 

+  1  ,  J+l  ,K.)  /DELY(  J)>   -2 


U5 


r 

1( 

G 

G 
G 
H 
U 


2d  =  V 

2C=v 

2=5e 

1A  =  V 
ib=V 
1C=V 

lO-r 

1=DF 
=  <E2 
=  (w( 

4.0E 

1  =  U( 

2  =  U( 

3  =  U( 
=  2.0 
=  (U< 
TEE  = 
c  ITE 


U+l  ,  JtKJMl  .0E0/DELY(  J)'   >  2  -1 .0  50/  D  ELY  ( J  )    >2) 
(  1+1, J+1,K)/RELY< J)v*  2 
F2G 

LY(  J)     DELY(  JP   (r2£+EzE-E2C)/(DELY( J)+DELY(J  )  ) 
Cl-lt  J+lfK-)/DELY(J)»r*2 

(1-1  ,  J,K).'  (1  .OE0/DELY(  J)     :  2-1  .0  EG/ D  EL  Y  ( J  )       2) 
(  I- it J+lf K) /D2LY( J)*' 2 
EIC 

LY(J)DELY(J)- (E1A+EIB-E1C)/ (DELY ( J  )+DELY(J )  J 
-rl)/(2.0^u-  CELX) 

I+ltJ,K)-W{I+l,JfK-l)-W(I-l,J,K)+w(I-l,J,K-l))/ 
0.  CELX"  DELZJ 
If J+1,K)/DELY(J) 

I  »  J,K)-  (1.0EO/DELYU  ) +1  .0  EO/ DEL  Y(  J  )  ) 
I  ,  J+l.K)  /CELY(J) 

E0'J  (G1-G2+G3)/(UELY( J  J+DELYU)  ) 
I  f  J,  K)-2.0E0*  U(I  ,  J,K)+U(  I,  J,K-1)  )/DELZ-3-  2 
-( A+B+C+D+REI" (E+F-G-H  )  ) 
(19,903)    UTEE 


CUCUL 
£=(U(I 

ai*v(  i 

B2=\M I 
B3=V< I 

B3=-E3 
B^DELY 
C  =  (W(I 
D=  (PHI 
E=(U( I 
1(  4.0FG 
F2A=V( 
C2E=V( 
F2C=V( 
F2C=-F 
F?=CEL 
F1A=V( 
C1E=V( 
F1C  =  V( 
F1C=-F 
FI=DEL 
FMF2- 
G=(W(  1 
HI =W( I 
H2=W( I 
Hj=W (  I 
H=2.UE 


A 
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+  1,  J 
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+G+C+R-I-(E+F-G-H) ) 
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ILCUL^TICW    OF     UT    FROM    3D    Ml 

ALCULATIGN    GN    A    CORNFR 

+1    AND     K-l    TERMS    HAVE    BEEN 

r    26     1-2, LEI 
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=  1 

*<U(  I 

1=U(I 

2=L)(  I 

5  =  U(  I 

=  DELY 

-cue  I 

=  (PHI 
2A  =  V( 
2d  =  V( 
2C  =  V( 
2=D£L 
1A-V( 

ie=v( 
1C-VC 

1  =  DEL 
*(E2- 
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4. OEO 
1  =  U(  I 
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3=U(I 
=  2.0E 
=  <U(I 
TEE=- 
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ME  N.TUM    EUU^TI  CM 
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,  J,K 

,  J-l 
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1+1, 
1+1, 
Y(  J) 
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I-lt 
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,J,K 
fJtK 
»  J-l 
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)/DE 
)  •  (  1 
,K)/ 
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K)/D 

Y(J- 

EO*D 
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LZ) 

LY(  J 

.OEO 

DELY 
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);  (B1+B2-B3 )/  ( DELY( J )+DFLY(J-l ) ) 

J,K+l)-U(I,J,K)nW(I,J,K))/(2.0E0lDELZ) 

I(  I-1,J,K )  )/( 2.0E0*DELX) 

Y(d»**2 

0E0/DELY(  J-D-   "  2-1.0E0/DELY(  J)'  -2) 

ELY(  J-1K-  2 

1)  '   (E2A+E2B-E2C)/(DELY( J)+DELY( J-i )) 

Y  ( J  ) "  :  2 

OEO /DELY  (J-l  )*■  2-1.  OEO/  DELY  (J  1**2) 

ELY(  J-lh  <2 

1 U  (E1A+E 16- E 1C)/(DELY( J)+DELY( J-l)) 

ELX) 

I  +  l,J,K)-WlI-i,J,IO-l)+w(I-l,J,K))/ 
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=  (LHI  +  l,J,K)*W(I  +  l,J,K)-U(I-l,J,K)W(I-l,J,K))/(2.QE0iDELX) 

1=V(  I,  J,«r.  Wtl,  J,K)/DELY(  J  Y**2 
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3  =  \A(  I,  J-1,K)*W(  I  ,  J-1,K)/DELY(  J-1)*'1  2 

=  DELY(  J-D'  DELYt  J  )*  (  B1+B2-B3  )/  (  CELYU  )+DELY(  J-l)  ) 

=  (W(  I  ,  J,K  +  1)*     2-*(  I,  J,K)  ;  ''  2)/(2.0E0*DELZ) 

=  (PhI  (  I,  J,K+1)-PHI  (  I,  J-,K  )  )/(2.0fcOfDELZ  ) 

MUd  +  1,  JtK  +  l)-U(  1+1  ,  J,K)-U<  1-1,  J,J<  +  1)  +U(  1-1,  J,  K)  )/ 

4. OEO'  DELX:  DELZ  ) 

2A=V{ I,J,K+1 )/DELY(J)**2 
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JfK)  (  1 
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03)     UT= 


U(I-l,Jf  K)  -2)  /(2.0E0    DELX) 

,K)/(DELY( J)     -  2) 

»  K)>  (1  .0E0/DELY(J-1)     •  2-  1  .GEC7D  EL  Y(  J  )  -'     2) 

, J-1,KI/D=LY< J-l)' l2 

1  )     (B1+B2-B3)/ ( DfcLYl J)+DELY< J-l) i 

,K)-U(  If  J,K-1  )=i  W(I,  J.K-1)  )/(  2.0  E0-- DELZ) 

HI(  1-1, J,K)  )/(  2.  GEO4  T:LX) 

LY(J)*  '2 

.  OEO/LELY(  J-l)*  > 2-1 .GEO/ C ELY (J  )"-  2  ) 

DELYC J-l)** 2 

-1)-:  (E2A+E2B-E2C)/(  DEL  Y(  J  )+DELY(  J-l)) 

L  Y  (  J  )  *  *  2 

.0E0/DELY(J-l)*r2-1.0E0/DELY(J)**2) 

DELY( J-l  )    *2 

-1)     (E1A+E  18-E 1C)/(DELY( J)+DELY( J-l) ) 

DELX  ) 

+1, Jf K-1)-W(I-1 f J,K)+W ( I-lf Jf K-l ) )/ 

J) 

0/DFLY(J-l)  +  1.0E0/DELY( J)  ) 

Y(J-l) 
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A=(U(I+ 
31=V(I , 

62=v(  i, 

B3  =  V(  It 
8=TELY( 
C=(W(  I, 
D  s  (  PH I  ( 
E  =  (U(  1  + 
1(4.0=0 
F2/i-V(I 
F2b="V(  I 
F2C  =  V(  I 
F2=0ELY 
F1A  =  V(  I 
FiB=V(I 
F1C=V(  1 
F1=0ELY 
F=(F2-F 
G=(W(  1  + 
H1=W(  I, 
H2=W(I , 
H3=W(  I  , 
H=2  .GEO 
WTEE=-( 
W^ITE(-2 
?   COM  I  Nli 


It  JtK 

J  t  K  M 

J,K) 

J-l  »K 

J-l^ 

J,KK 

It  J,K 

It  JtK 

CFLX 

tJtK) 

.JtK) 

tJ-lt 

(  J)  >? 

iJfK- 

.JiK- 

tJ-lt 
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It  J»K 
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J-1,K 
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W(I t JtK 
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DELY( J) 
'  2- W  (  I , 
)-PHI ( 
)-U(I+ 
DELZ) 
/CELY( 
Ml.OE 
K  )/DELY 
ELY( J-l 
1 )  /  D  E  L  Y 
DM1.0 
K-l) /DE 
ELY( J-l 
.OEO*DE 
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DELY(J  ) 
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)/DELY( 
H2+H3  )/ 
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,J,K)-U(I-1,J,K)W(I-1,J,K) )/(2.0E04  DELX) 

)/DELY(  J)*'»2 

)     ( l.OEO/OE LYf J-l) »     2-1.0E0/DELYU)* -2 ) 

-1,K)/DELY( J-l)>-  2 

1  CB1+B2-B3)/  (DELY(J)-rCELY(J-l  )  ) 

J,K-1)  *■.  2)/(  2.0E0±DFLZ) 

,J,K-1 )  )/(2.0F0-  DELZ  ) 

t  J,  K-l) -1/(1 -1  t  JtK)+uri-l  tJtK-1)  )/ 

)     '  2 

/DELY(J-l)        2-1.0E0/DELY(J)-;<2) 

(J-  l)*-2 

)>  (F2A+F2B-F2C)/(DELY(  J  ) +OELY  (  J-l  )  ) 

(  J  ) ;      2 

/DELY(J-l)^  2-1.0/ DEL Y( J  )*' 2 ) 

LY(  J-l)  '•  ^2 

)•  (P1A  +  F1B-F1C)/(DELY(  J)+DELY(  J-D) 

LZ) 

,J,K))  /(2.0E0^DELX) 

DELY(J-1)+1.0E0/DELY ( J)  ) 

J-l) 

(DFLY(J )+DELY(J-l) ) 

(E+F-G-H) ) 


2 
C 
C  CALCULATION  AT  ENDS  OF  THE  REGION 


2£ 


DO  28  J=l f ME 
GO  28  K=1,NF 
1  =  1 

UTEE=U. OEO 
WRITE(  1*7,903) 
W^EE=O.CEO 
WRITE(20,9G3) 
CONTINUE 


UTEE 
WTEE 


MC 


J=l, 

K=1,NE 


DH  2^ 
DO  29 
I=LE 

U"!"EE=O.OEO 
WFITE(19,903) 
wTEE=0. CEO 
WRITF(20,9C3) 
CONTINUE 

IF( ITER .EO.O) 
IF (  ITER  .EQ.l ) 


UTFE 

WTEE 

DELT=DELT2 
DELT=CELT1 


h9 


^0? 


IF  C  IT  EF  .cQ.O)    GG    TC    402 

IfMICALC.EO.O)     S4j    TG    402 

REWlNC     16 

R  E  W I N  C     18 

RFA0C16J     U 

R  C  A  D  (  1  b  )     W 

CCMINUF 

REHINC     19 

REWINC    20 


410 


5  FAD    VALUES    GF    UTCE    ANC    WTEE    AND    COMPUTE    NEW    U    AND    W 

DC    410    1=2, LEI 

Dr>    410    J=2,ME1 

DC    410    K=2,NE1 

kEAu(19,903)     UTc" 

U(  If JtK)  =  U( I, J,K)+UTEE'  DELT 

READ(20,904)    WTEE 

h(  I  ,J,K )=W( I  , J,K)+WTEE;  DELT 

CCNT  IMJE 


411 


DG  411  1=2,  LEI 
DC  411  J  =  2tv<El 
K=l 

P::AD(  19,903)     UTEE 
U(I,J,K)=UII,J,K  )+UTEE»  DELT 
SEAr)(20,904)     WTEE 
W(l,J,K)  =  *(I,J,K)  +  l*TEE*DELT 
CXNT  INUE 
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EC     412    1=2, LEI 

OG    412    J=2,ME1 

K*NE 

F<EAD(19,903)     UTE£ 

Ul  I,J,K  )=U(1,  J,K)+UTFE- 

kEAD(2C,904)     wTEE 

W  (  I, J,K  )  =  W( I,  J,< J+toTEE 

CENT INUE 


DELT 

DELT 


413 


DC    413 
Cf    413 

J  =  l 

^EAD(19,903)     UTEE 

U(  I  ,J,K)  =U(I  ,  J,K)+UTEE:rDELT 

WT^F 

J,K  l+WTEE*  DELT 


1=2, LEI 

K=2,^E1 


9  04) 


CEAD( 20, vun* 
rt( It Jv K  )=W(  I, 
CCNTIMUE 
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414 


415 


416 


417 


418 


DO    414    1=2, LEI 

CO    414    K=2,NE1 

J  =  ME 

PEAD(19,903)     UTEE 

U(  I,J,K  )=U(  I,J,K  )+UTEF^ 

RF£D(20,904J     WTEE 

W(  I,J,K  )  =  W(  I,JfK)+WTEE: 

CONTINUE 


DELT 
CELT 


DO    415     I  =  2,LE 

J=l 

K  =  l 

kEAD(19,903) 

U(  I,J,K)=U( I, 

READ(20,904) 

W(I,J,K)=W(I, 

CONTINUE 


UTE  = 

J,K)+UT  EE*OELT 

WTEE 
J,K)+WTEE-DELT 


1=2, LEI 


DC    416 
J  =  l 
K  =  NE 

READ(19,903) 
U(I  ,J,K)=U(I 
READ(20,904) 
W( I, J,K )=W(I 
CONTINUE 


UTEE 

J,K)+UTEE~ 
WTEE 
J,K)+WTEF" 


DELT 
DELT 


CO    417     1=2, LEI 
J  =  K'E 
K=l 

PEAD(19,903) 
U( I ,J,K)=U( I 
PEAD(20,904) 
W(  I,J,K  )  =  W(I 
CONTINUE 


UTEE 

J,K)  +  UTEE" 
WTEE 
J,K)+WT  EE 


DELT 
DELT 


DC    418     1=2, LEI 

J=ME 

K  =  NE 

READ(19,903) 

Ul I,J,K)=U(I 

READ(20,904) 

W(  I  ,  J,K) =W(I , J,K)+WTEE 

CONTINUE 


UTtE 
J.KJ+UTEE* 

WTEE 


DELT 
DELT 


DC  419 
DO  419 
1=1 


J=1,ME 
K=1,NE 
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c 

c 
c 


U (I, J,«  )  =  1.0E0 
M I ,J,K) =0.0E0 
4  19    CONTINUE 

Dr    420    J=lfME 

DO    420    K=1,NE 

U(l.EfJf<)  =  3.0E0-U(LEl,JfK)-3.0£O'U(LE2tJ,K)+J(LE3,J,K) 

W(L  =  ,J  ,K)=3.0E0*W(  LEI  t  J~,  KJ-3.3E0    w(LE2t  J,K)+w  (t.^3,  J,K) 

420  CONTINUE 

DC    421     I=LP1  ,LP2 
CC    421    K=:MPlfM°2 
J  =  l 
U( 1 ,JtK] =0.0E0 

w( if j,k  )  =  o.oeo 

421  CONTINUE 

CALCULATE    V    WITHIN    REGION    FROM    3D    CONTINUITY 

DC    30    I =2, LEI 

lr     30    J=1,ME1 

CO    30    K=2fNEl 

UCX2=(  U(I+lf J,K-1 ) +U(  1+1, J  +  lfK-1  )     +U(I  +  l,J,K  +  l)+U(I  +  l,J-t-l,K  +  l))/ 
14.0E0 

uEXl=(U(I-lfJ,K-l)+U(I-lfJ  +  l,K-l)+U(I-l,J,K+l)+U(I-l,J  +  l,K+in/ 
14. CEO 

JEX=(UEX2-UEXl)/(  2.0E0-DELX) 

v.Z2  =  (  W(  I+lt  J  fK  +  1  )+W(I+lt  J  +  ltK+l)+Wf  I-li  Jf  K+l)+W(  1-1,  J+lfK-f-UJ  / 
Ih.OEO 

«Z1=(W( I+l,Jf K-l)+W( I+1,J+1,K-1)+W(I-1, J,K-l)  +  w( I -If J+l ,K-1) ) / 
14.0E0 

\nlE-(  WZ2-WZl)/(  2.  OEO-CELZ) 

»M  1  f  J+lf  K)=V(  I,  J,K)-DELY(J  )     (UEX  +  WZE) 
30    CCNTINUE 

CALCULATE    V    ON    A    SIDE    FROM    3-D    CONTINUITY    EQUATION 

K-l    TECMS    WVE    SEEN    ADJUSTED 

DO    31    I=2fLEl 

DC    31    J  =  1 , Mc 1 

K  =  l 

UEX2=(U(I+lfJ,K)+U(I+l,J+lfK)     +  U(I  +  l,J,K+l)+U(I  +  l,J  +  lft<  +  l))/ 
14.0  EC 

UEXl=(U(I-l,JfK)+U(I-lfJ+lfK)+U(I-lfJ,K+l)+U(I-l,J+lf<+l))/ 
14.0E0 

JEX=( UEX2-UEXl)/(2 .QEO    DELX ) 

viL  ?  =  (  '*(  I  +  l,Jf  K+1)  +  W(I  +  1,J+1  ,K  +  1)+W(  1-1  ,  J,K+1)+W(I-1,  J+1,K+1)  )/ 
14.0E0 

WZ1  =  (  w(  1  +  1  ,  JfK)  +-W(  1+1  t  J  +  l  flO+rt  (  I -If  J,K)  +  W  (  1-1,  J  +  l,  K  )  )/ 
1  4  .  C  E  0 
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WZE=(  WZ2-WZ1) /( 2.0E0- DELZ) 
V(  I,  J  +  l,k)=V( I, J,K)-DELY(J  )  MUEX+WZE) 
31    CONTINUE 

r 

C  CALCULATE    V    CM    A    SICE    FRDM    3-D    CONTINUITY    EQUATION 

C  K+l    TERMS    AHVE    BEEN    ADJUSTED 

DO  3  2  1=2, LEI 
CO  32  J=1,ME1 
K  =  NE 

UEX2=(U(  1  +  1,  J,K-1)+U<  1  +  1,  J+lt  K-l)     +U(  1  +  1,  J,  K) +11(1  +  1  ,J+1,K)  )  / 
14. CEO 

UFX1={  U( I-i, J,K-1)+U( 1-1, J+1,K-1 )+U(I-l  , J,K)+U( 1-1, J+l, K) )/ 
14.0EO 
UEX  =  (UEX2-UEX1 )/(2 .0 E0>  DEL X  I 

WZ2=U( I+1,J,K)+W{ 1+1 , J+1,K)+W( 1-1, J,K) +W(I~1, J+1,K) )/ 
14.0EC 

WZ1  =  (  W(  1  +  1 ,J,K-1)+W(I+1 ,J+1,K~1 )+W(I-l , J, K-l  )+w(  1-1, J+1,K-1)  )/ 
14.0E0 
WZE  =  (  WZ2-WZ1  )/(2.0E0    DELZ) 
V(  I,J+1,K)=V( I , J,K)-DELY( J)>  (UEX+WZE) 
3  2    CONTINUE 
C 

C  SET    V    EQUAL    TO    ZERO    IN    THE    PL£NE    GF    THE     PLATE 

DC    34     1=1, LE 
DC    34    K=1,NE 
V(  I,  1,K)  =  0.0E0 
34    CONTINUE 
C 

C  SET    UPSTREAM    V    EQUAL    TO    ZERG 

DC    35    J=1,MF 
DO    35    K=1,NE 
V(  1,J,K  )  =  0.0E0 
3  5    CONTINUE 
C 

C  EXTRAPGLATE    DOWNSTREAM    V    FIELD 

DO    36     J=l ,  ME 
DC     36    K=1,NE 

V(LE,J,K)  =  3.0E0*V(LE1, J,K)-3.0E0<  V(LE2, J,K)+V(LE3, J  ,  K) 
36  CONTINUE 
C 

WRITE (6, 8 12) 
£12  <=CFMAT(  1H0,20X,'  INTO  CALCULATION  POISSON') 
200  CONTINUE 
C      CALCULATE  PRESSURE  FIELD  FROM  PGISSONS  EQUATION 

EK.R*0.  OEO 
C 

C  CALCULATION    OF     PRESSURF    FROM    3D    POISSON'S    EQUATION 

C  CALCULATION    WITHIN    REGION 
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c 
c 
c 
c 


37 


DC  37 
C~  37 
CC  37 
A  1=(P 
A2=(P 
A  3  =2. 
KOELY 
A=A1  + 
3  =  (U( 
C1=V( 
C2=(V 
C3  =  V( 
C=2.0 
D  =  (W( 
=  2A=U 
F2E=U 
E2C=U 
E2=DE 
E1A=J 
E1B  =  U 
E1C=U 
E1  =  DE 
E  =  (E2 
F2A=V 
F2B  =  V 
^2C=V 
F2  =  DE 
F1A  =  V 
F1R=V 
F1C  =  V 
F1=DE 
F=(F2 
G2=U( 
G1=U< 
G=(G2 
H  =  (2. 
APHI  = 
ft  E  S  I D 

FCC  =M 

PHKI 
CCNTI 

CALCJ 
CALCU 
K-l  T 
D:  38 
DC  38 
K=l 


1=2, 
J  =  2, 
K=2, 
HI  (  1  + 
HI  (  I, 
OEO*  ( 
(J  )  +  D 
A2+A3 
I  +  lt  J 
It  J  +  l 
(I, J, 
I  t  J-l 
EO-MC 
I  tJ,K 
(I  +  lt 
(I+lt 
(I+l, 
L  Y  (  J  ) 
(1-1, 
(  1-1, 
(  I-lt 
LY(  J) 
-El)/ 
(  ItJ  + 
(I  tJt 
(  I  ,J- 
LY  (J) 
(I  ,J+ 
(  I  tJt 
(  I  ,J- 
LY(  J) 
-Fl)/ 
I  +  lt  J 
I-lt  J 
-Gl)/ 
OEO/D 
(  l.OE 
=  ABS  ( 
AXKE 
tJfK) 
NUE 


LEI 
ME1 
NE1 

1,  J,K)+PHI ( 
J,K+1)+PHI  ( 
PHI (I t J+ltK 
ELY(J-l)  ) 

,K)**2-2 .OE 
,K  )■:•■■  2/DELY 
K.  )*  "2)  *(  1.0 
,K)*"  "2/DELY 
1-C2+C3)/(D 
tl)*>2-2.0E 
J+ltK)'  VU  + 
J,  K) -\f(  I  +  l, 
J-lfK)  •  V(  1  + 
" DELY( J-l) 
J  +  ltK  )*V(  I- 
JtK)*  Vll-l , 
J-1,K)*V(  I- 
-DELY( J-l) 
(  2.0E0*DELX 
1,K  +  1)'*W(I, 
K+1)*W( I t Jt 
ltK+1)  vw(  I  , 
■'DELY(J-1)> 
ltK-lJ*W( It 
K-l)-  W.(  I  ,J, 
I,K-1)*W(  It 
»  DELY(  J-l)* 
(2.0EO-DELZ 
,K+1)*W( I+l 
,K+1)<  W(I-i 
(4.0EO-*DELX 
EL X* ;  2+2  .OE 
O-GMEGA )-PH 
APHI-PHI ( I, 
DR,RESID) 
=  APHI 


I-lt Jt 

It JtK- 

)/DELY 


0*U(I, 
(J) 

EO/DEL 
(J-l) 
ELY(  J) 

0  « W  (  I  t 

1  »J+i  t 
J  ,  K  )  v  ( 
1 tJ-1 t 
(E2A+E 
ltJ+lt 
JtK)  M 
ltJ-lt 
(ElA+E 
) 

J  +  1.K  + 
K  +  l)  ( 
J-ltK+ 
(F2A+F 
J+ltK- 
K-l)  --( 
J-1,K- 
JF1A+F 

t JtK+1 
t J, K+l 
>  DELZ) 
O/DELX 
KI  ,Jt 
JtK)  ) 


K))  /DELX'  A-2 
1))/DELZ**2 

(J)+PHI U,J-1,K)/DELY( J-l)  )/ 


J , K)**2-rU  ( 1-1 1  J  t  K)**-2 ) /  DELX*    2 

Y( J-l)+i.O/DELV( J) ) 

+  CE  L  Y (  J—  1 )  ) 

J,K)**2+W<  It  J,K-l)*-2)/DELZ--*2 

K)/DELY(J)*     Z 

1.0EO/DELY(  J-l)>     2-1.  OE  O/DE  LY  (  J)  *  4  2) 

K-)/D  ELY  (J-l  1**2 

26-E2C)/(DELY( J)+DELY( J-l)) 

K  )  /  D  E  L  Y  (  J  )  *  *  2 

l.OEO/DELYU-1  )**2-1.0E0/DELY(J  )       2) 

K)  /DELY(  J-l)  v»  2 

1B-E  1C)/(DELY(J  )+DELY(  J-l)) 

1)  /DELY(  J)^2 

1 .0  EG/ DELY ( J- l)^'2-1.0EO/ DEL Y( J  )    -  2) 

1)/DELY( J-l)  ;*2 

23-F2C)/(DELY( J ) +DEL Y( J-l )  ) 

1)/DELY(J)**2 

1.0EO/DELY(  J-l)*»-2-1.0F0/DELY(  J)--2  ) 

l)/DELY(J-l)**  2 

18-F1C)/(DELY(J)+DELY(J-1  )) 

)+U(  I-lt  J  t  K-l  )-W(I-l,  J,  K-l) 
)+U(  I  +  lt  JtK-Ir)':  U  (I+l,  J,  K-l) 

*  ■  -2+2.0E0/(CELY(J)*DELY(J-l)  )  ) 

K)+CMEGAr (A+B+C+D+2.0E0- (E+F+G) )/H 


LAT  ION  OF 

LATICN  ON 

ER.MS  HAVE 

1=2, LEI 

J=2,ME1 


PRESSU 
A  S  IDE 
BEEN  A 


RE    FRO 
DJUSTE 


M    3D    PCISSONi'  S    EQUATION 

D 


5h 


c 
c 

c 
c 
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A1=(PHI ( I+1,J,K) 
A2=(PHI ( It  J,K>1) 
43=2. OEO* '(PHI  (It 
1(DELY(J  )+DELY(  J- 
A=AH-A2  +  A3 
b  =  (U(  1  +  1,  J,K)*t2 
C1=V( I f J+I,K)**2 
C2=(V( I ,J,K)**2) 
C3=V( I, J'1,K1**2 
C=2.0E0MC1-C2+C 
D=(  W(  I  ,J,K+1)^2 
E2A=J( 1+1, J+1,K) 
E2B=U(I+1, J,K)*V 
E2C=U( 1+1, J-1,K) 
E2  =  DELY(J)*"DELY( 
ElA=U(I-lf J+l fK) 
E16=U(  I-1,,J,K)*V 
E1C=U( 1-1, J-1,K) 
£1  =  DELY(  J)OELY( 
E=(E2-E1)/(2.0E0 
F2A=V(I  ,J+1,K+1) 
F2B  =  V1 I  ,J,K+I)*W 
F2C=V( I ,J-1,K  +  1) 
F2=DELY(J)*DELY( 
F1A=V( I,J+1,K)*W 
F1B=V( I ,  J,K)*W( I 
F1C  =  V(  I  ,J-1»K)*W 
F1=DELY(JHDELY( 
F=(F2-F1)/  (2. GEO 
G2=U(  I+1,J,KH)  * 
G1=U(I-1,  J,K+1I* 
G=(G2-G1)/(4.0E0 
H=(2.0E0/DELX**2 
APHI=(1 .OEO-OMEG 
RESIQ=ABS(APHI-P 
EPP  =  MAXKERRf  RES 
PHIM  ,  J,K)=APHi 
CONTINUE 


+  PHK  1-1,  J?K))/DELX**2 

+PHI (I ,J,K) ) /DELZ**2 

J+1,K)/DELY( J)+PHI(I , J-1,K) /DELY( J-l))/ 

1)  ) 


-2.0E 
/DELY 
*(1.0 
/DELY 
3)/(D 
-2.0E 
*V(I  + 
(  1+1, 
*V(I  + 
J-l)* 
*V(I- 
(1-1, 
*V(I- 

J-ll* 

*DELX 
*  W  (  I , 
(I  tJf 
*W  (  I , 
J-l  J* 
(  I  ,J+ 

t  Jfio 
(i  ,j- 

J-l)* 
i-DELZ 
WU  +  1 
W(  1-1 
*DELX 
+  2.0E 
A)*PH 
HI  (I  , 
ID) 


0*  U(  I 

(J) 

EO/DE 

(J-l) 

ELY(J 

0*W(I 

1,J+1 

J,K)* 

1,J-1 

(E2A  + 

If  J+l 

J,K)* 

If  J-l 

(E1A+ 

) 

J  +  lf  K 

K+l)* 

J-tfK 

(F2A  + 

1,K)/ 

*  (1.0 

1,K)/ 

(F1A  + 

) 

f  JfK  + 

f  JfK  + 

*DELZ 

O/DEL 

I  (If  J 

Jf4<)  ) 


,  J,KM*2-HJ(I- 
LY(J-1X+1.0/D 


)+DEL 
f  JfK) 
f  K)/D 
(l.OE 
,K)/D 
E2B-E 
f  K)/D 
(  l.OE 
,K)/D 
E1B-E 


Y(J-l)  ) 
**2+W(If 
ELY! J)** 

0/DELY(J 
ELY( J-l) 
2C)/(DEL 
ELY(J)** 
0/DELY( J 
ELY(J-l) 
1C)/(DEL 


1  , J,K)**2)/DELX**2 
ELY(J)  ) 


J,K)**2)/DELZ**2 

2 

-1)**2-1.  OEO/DELYC  J  )**2) 

**f  2 

Y( J)+DELY(J-1)) 

2 

-1)**2-1.0E0/DELY(J)**2) 

**  2 

Y(J)+DELY(J-D) 


+  1  )/DELY(J)** 
(1.0E0/DELY(  J 
+1)/DELY( J-l) 

F2B-F2C)/(DEL 
DELY( J)**2 
E0/DtLY{J-l)* 
DELY ( J-l )**2 
F  1B-F1C)/(DEL 

1)+U(I-1,J,K) 
1)+U( I+lf J,K) 
) 

X*A2+2.0E0/(D 
,K)+GMEGA*(A+ 


-1  )**2-l. OEO/DELYC J)** 2) 
Y(J)+DELY(J-D) 
*2-1.0E0/DELY(  J)#*2) 
Y( J)+DELY( J-l) ) 

*W( I-1,J,K) 

*W(I+1,J,K) 

ELY(  J)*DELY(J~-1)  )  ) 

B+C  +  D+2.0E0*(E+F  +  G)  )  /H 


CALCULATION    QF    P 
CALCULATION    ON    A 
K+l    TERMS    HAVE    B 
DO    39    I =2, LEI 
DO    39    J=2,ME1 
K  =  NE 

A1=(PHI (I+l,JfK) 
A2=(PHI  (  I,.JtK)+P 
A3=2.0E0*(PHI (I , 


RESSURE    FROM    3D    POISSON'S     EQUATION 

SIDE 
EEN    ADJUSTED 


+  PHI  (1-1, J,K))/DELX**2 
HI( I, J,K-1) )/DELZ**2 
J+1,K)/DELY(J)+PHI (I,J-lfK)/DELY( J-l) )/ 
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c 
c 
c 
c 


39 


1  U£LY( 

A=A1+A 

B=(IMI 

C1=V(I 

C2=(V( 

C3=V(I 

C=2.0E 

D=(W(I 

E2A=J< 

E2B=U( 

E2C=U( 

E2  =  DEL 

E1A=U( 

E1B=J( 

E1C=U( 

rl=DEu 

E=(E2- 

F2A=V( 

F2b=Vl 

F2C=V( 

F2=DbL 

FIA=VI 

F1B*VI 

FiC=V( 

P1=DEL 

F=(F2- 

G2=U( I 

G1=U(I 

G=IG2- 

H=(2.0 

APHI=( 

RESID= 

ERR=hA 

PHI  (I  , 

CCNTIN 

C^LCUL 

CALCUL 
J+l    TE 
DC    40 
DC     40 
J  =  ME 
A1=(PH 
A2=(PH 
A3=2.0 
1( DELYI 
A=A1+A 
BMUU 


J  )+DE 

^+A3 

+  lf  J, 

.J+l, 

I  f  JfK 

t  J- It 

0MC1 

t  JtK) 

1+1,  J 

I+lf 

I+lf 

Y(  J) 

I-lt 

I-lt 

I-lt 

Y  (  J  )  * 
El)/( 
I  ,J  +  1 
I  t  JtK 
1  ,  J-l 
Y(  J)" 
I  ,  J+l 
I  .JfK 
I  ,J-1 

Y  (J)* 
Fi)/{ 
+  lt  Jt 
-It  J, 
Gl)/< 
EO/DE 
l.OEO 
ABS(  A 
XKER 
J  ,  K )  = 
UE 


LY( J-l  )  ) 


KH*2- 
K)**2/ 
)**2)"* 
K)**?/ 

-C2+C3 

**2-2  • 

+  1, K  )  r 
tK)1-V( 
-ltK)* 
DELY( J 
+1,K)* 

t  k  i  *  \r< 

-1,K)* 
DELY(  J 

2  .OEO^ 
,  K I  ♦  W  ( 

)*  Wl  I  , 
tK) -w ( 

3ELYI J 
t  K- 1 ) * 
-1)» W( 
,K-1)  f- 
DELYIJ 
2.0EC* 
K)iW(  I 
K)*M  I 
4.GE0J- 
LX**2+ 
-OMEGA 
PHI-PH 
Rt RESI 
APHI 


2. 
DE 
(  1 
DE 
)/ 
OE 
VI 
1  + 
VI 

-1 

V( 

I- 

VI 

-1 

CE 
I  » 
Jt 
I  t 
-i 
W| 

If 

W( 

-1 

DE 
+  1 
-1 
DE 
2. 
)* 
I  ( 
Dl 


OEO*U( 
LY  (J) 
.OEO/D 
LY(J-1 
(CELY( 
0*  W 1 1  f 
I+1,J+ 
If  JtK) 
I+lf J- 
)*  (E2A 
I-1,J+ 
It  JtK) 
I-1,J- 
)i  (E1A 
LX  ) 

J+l ,K) 
M* (  1. 
J-ltK) 
)" (F2A 
It  J  +  lt 
JtK-i) 
I tJ-lt 
)*  (F1A 
LZJ 

fjf K)+ 
,  JtK)  + 
LX*GEL 
0E0/D5 
PHI  (I, 
I ,JtK) 


I  t  J=VK)'*i2+U(  I-1,J,K)**2)  /DELX*  *2 

EIY(J-1)+1.0/DELY( J) ) 

) 

J)+DELY(J-1)) 

JfK)**2+hCI ,J,K-1)**2)/DELZ^2 

If K)/DELY|JJ**2 

Ml.OEG/CELYIJ-1  K^-l  .OEO/DtLYU  )**2) 

If  KJ/DELYI  J-l)*-*  2 

+  E2B-E2C)/(DELY( J )+DELY(  J-l)) 

1,K)/DELYIJ)**2 

•I l.OEO/DELYt J- 1)^*2-1. 0E0/G5LY( J)^*2) 

1,K)/DELY(J-1)*~2 

+E1S-E ICi/IDELYI J)+DELY( J-l ) ) 

/DELY( J)**2 

OEO/DELY(  J-1)M  2-1.0E0/DELY(J)**2) 

/DE4_Y(  J-l)-*v2 

+F26-F2C)/(DELY(J)+DELY( J-I) ) 

K-D/DELYC  J)**2 

* (1 .OEO/DELYlJ-1  )*42-1.0E0/D£LY(J )**2) 

K-1)/DELYU-1)**2 

+F1B-F1C)/(DELY( J)+DELY( J-l)) 

U( I-lf JfK-l)*WlI-lfJfK-l) 

U(I+l,J,K-l)*W(I+lfJ,K-l) 

Z) 

LX**2+2.0E0/(  DELYt  J)*>DELY(  J-l)  )  ) 

J,K)+OMEGA*(A+B+C+D+2.0EO*( E+F+G) )/H 

) 


ATI CN    OF 
ATIOi\    ON 
RMS     HAVE 
I  =2,  LEI 
K=2,NE1 


PRESSURE    FROM    3D    POISSON'S     EQUATION 

TO  P 

3EEN    ADJUSTED 


I  I I  +  lf J,K)+PHII  1-1, J,KJ) /CELX**2 

I I  If JfK  +  l)+PHI tit JtK-l))/DELZ**2 

EGM  PHI  (I,  J,K)/DELY(J  )+PHI(I,J-l,K  )/DELY<  J-l)  )/ 

J)+DELY(J-1) ) 

2+A3 

+  1,  J,K)  M2-2  .GEO*  Ut  I  ,  J  ,  K  )  -.'J/2  +U  ( 1-1 ,  J  ,  K  )**2  )  /  DELX*4  2 
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C1=V(I , J,K>**2/DELY(J) 

C2=(V( I,J,K)**2)*< 1.0E0/DELY( J- 1 ) +1 .O/DE LY ( J) ) 

C3=V(I,  J-ltK)**-2/DELYU-l) 

C=2.0E0MC1-C2  +  C3)/(DELY(  Jl+DELY (J-l)) 

D=(Wl  I  ,  Jt  K  +  1K*2-2.0E0*W(I  ,J,K)**2+W(I,J,K-1)**2)  /DELZ**2 

E2A=U( 1+1, J,K)*V< I+l» J,K)/DELY( J)**2 

E2B=U(  1+1,  J»K)*V(I-H,  J,K)  M1.0E0/DELY{J-1)**2-1.0E0/DELY(J  )**2) 

E2C  =  U(  1+1,  J-1,K)*V(  I  +  1,J-1,K)  /DELY(J-1)**2 

E2=DELY(J)*DELY< J-l)*  (  E2 A+E2B-E2 C )/ ( DEL Y < J  )+DELY( J-l)) 

E1A  =  U(  I-1,J,K)*V(I-1,J,K)/DELY( J)**2 

E1B=U(  1-1,  J,K)*V( 1-1, J,K)*{  1.0E0/DELY(  J-l)**2-1.0E0/DELY(J)-**2) 

E1C=U(I-1,J-1  ,K)*V(I-1  ,J-l,K)/DELY(J-l)**-2 

E1=DELY(J)*CELY< J- 1 )*  (  E 1A  +  E 1B-E 1C) /( DE LY ( J) +DELY ( J-l ) ) 

F=(E2-E1)/(2.0E0*DELX) 

F2A=V( I ,J,K+1)*W(I  ,J»K+1)/DELY( J)** 2 

F2B=V( I  ,J,K+1)  M(  I, J,K+1)*(  1.0E0/DELY(  J- 1 )** 2-1 . OEO/DE LY (  J)** 2) 

F2C  =  V(I  tJ-lf  K+l)*W(If  J-1,K+1)/GELY<J-1)**2 

F2=DELY( J)*DELY( J-l)*  (  F2A+F2B-F2  C)  /(  DEL  Y  (  J  )  +DELY  (  J-l )  ) 

FlA=V(If  JfK-l)*WUtJ,K-I)/DELY<  J)** 2 

F1B  =  V(I  ,J,K-1)*W(  I  »  J,K-1)*  (l.OEO/CELYCJ-1  )**2-1.0EO/DEL  Y(  J  )**2) 

F1C=V(I ,J-1»K-1)*W( I, J-lf K-1)/DELY( J-l)** 2 

F1=DELY(J)*DELY(J-1)*(F1A+F1B-F1C)/(DELY( J)+DELY( J-l) ) 

F=(F2-Fl)/(2.0EO*DELZ) 

G2  =  U(I  +  1,  J,K+1)  *W(  I-H,J,K+1)  +  U(I-1,J,K-1)*W(I-1,J,K-1) 

G1=U(  I-1«J,K+1)*W( I-lf JfK  +  l)+U(  I+lf JfK-1  )*W(I  +  lf Jf K-1J 

G  =  (G2-G1)/(4.0E0*DELX*DELZ) 

H=(2.0EG7DELX**2  +  2.0EO/DELX**2+2.0EO/(QELY( J)*DELY( J^l)  )  ) 

APHI  =  ( l.OEO-OMEGA)*  PHI ( I , J , K )+OMEGA* ( A+B+C+D+2.0E0* ( E+F+G) )/H 

RESID=ABS<  APHI-PHH  I,~J,K)  ) 

ERR=MAX1(ERR, RES  ID) 

PHKI  ,  J  ,K)=APHI 
40    CONTINUE 
C 

C  CALCULATION    OF    PRESSURE    FROM    3D    PCISSON1 S     EQUATION 

C  CALCULATION    DN    THE    BOTTOM 

C  J-l    TERMS    AHVE    BEEN    ADJUSTED 

DC    41    1=2, LEI 

DO    41    K=2,NE1 

J  =  l 

A1=(PH1(I  +  1,J,K)+PHI(I-1,J,K))  /DELX**2 

A2=(PHI  (I,JfK-H)+PHI(If  J,K-l))/DELZ+*2 

A  3=2.  GEO*  (PHI  (I  ,  J+l  ,K)  /DELY(  J)  +  PHI  H  ,J+1  ,K)/DELY(J  )  \J 
HDELVt  J  J+DELYU)  ) 

B  =  (U(I  +  1,J,K)*^-2.0E0*U(I  ,J,K)*+2+U(I-l,J,K)**2)/DELX**2 

C1=V(I, J+1,K)**2/DELY(J) 

C2=(V(I , J,K)**2)* (1.0EO/DELY(J)+1.0/DELY(J ) ) 

C3=V< I, J+1»K)*^2/DELY( J) 
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c 
c 
c 
c 


C  2=-C  3 

(>2.0F0*  (Cl- 
D=( W(I , J,K+1 
E2A=U( 1+1, J+ 
E2E=U(  1+1,  J, 
E2C=U( 1+1 ,J+ 
£2C=-E2C 
E2=DELY(J) *D 
Ei£  =  Ul I-1,J+ 
blB=Ul 1-1, J, 
E1C=U(I-1,  J+ 
E1C=-E1C 
Ei=DELY (J)*D 
£=(E2-El)/(2 
F2A=V( I ,J  +  1, 
F2B=V(I  ,J,K+ 
F2C=V( I , J+l, 
F2C=-F2C 
F2=DELY(J) *D 
F1A=V( I, J+l, 
F1B=V( I ,J,K- 
F1C  =  V(  I  ,J+1, 
F1C=-F1C 
F1=DELY(J)*-D 
F  =  (F2-Fl)/(2 
b2=U( 1+1, J,K 
G1=U(I-1, J,K 
G=(G2-Gl)/(4 
H=<2. OEO/DEL 
APHI=(  l.OEG- 
RESID=AbS(AP 
ERR=MAX1(ERR 
PHK  I  ,J  ,K)  =A 
41    CONTINUE 

CALCULATION 
CALCULATION 
J-l    AND    K-l 
DC    42    1=2, L5 
J  =  l 
K  =  l 

Al=tPHI  (1+1, 
A2=(PHI (I , J, 
A3=2.0EC*  (PH 
1  (CELYt  J  )+DEL 
A=Al+A2+A3 
B=(U(I+1, J,K 
C1=V( I, J+1,K 


ELY  (J)MB2A  +  E2B-E2C  )/(DELY(  J )+D 
l,*)?V(I-lfJ+l,K)/DELY{J)**2 
KMVU-1, J,K)*( 1.0fcO/DtLY( J)**2 
1,K)*  VCI-1  ,J  +  1  ,K)/DELY(J)»»-*-2 


1)*>2)/DELZ**2 
-l.O-EO/DELYl  J)**2) 

ELY(J)  ) 
-1.0E0/DELY(J)-»*2) 


ELY<  JM  (elA+ElB-ElC)/(DELY<  J  )+DELY(J)  ) 

.OEO*DELX) 

K  +  1)*WU,  J  +  1»K+1)/DELY(J)**  2 

1)*W( I , J,K+l)*(1.0EO/DELY(J  )**2-l  .OEO/  DELY(  J)**2) 

K+1)*W(I,J+1,K+1)/DELY< J)**2 


ELY( J)*(F2A+F2B-F2C1/ (DELY( J )+C 
K-1)*W(I,J+1,K-1)/DELY(J)**2 
l)*W(I,J,K-l)*(1.0EO/DELY(J)**2 
K-1)*W(I ,J+1,K-1)/DELY(J)**2 


ELY(J  )  ) 

-1. OEO/DEL Y( J)**2 


ELY( 
.  OFO 
+  1)* 
+1)* 

.OEO 
X**2 
CME3 
HI-P 
,RES 
PHI 


J)* (F1A+F1B-F1C)/(DELY(J )+DELY(J) ) 

^DELZ) 

W( 1  +  1, J,K+1)+U(  I- 1, J, K-l)*  V, 

W(I-1» J,K+1)+U( I+1,J,K-1)*W 

*DELX*DELZ) 

+2.0£0/DELX**2+2.0E0/(DELY( 

AH  PHI  (I , J,K)  +  OKEGA* (A+B  +  C  + 

HI(  I, J,K) ) 

ID) 


(I-1,J,K-1) 
(I  +  1,J,K-1) 

J}*D5LYUJ) ) 

C+2.0E0ME  +  F  +  G)  >/ 


H 


OF    PRESSURE    FROM    3D    POISSON'S    EQUATION 
ON    A    CORNER- 
TERMS    HAVE    BEEN    ADJUSTED 
1 

J,K)+PHI( 1-1, J,K)) /DELX**2 

K+D+PHK  rf  J.KJ  r/DELZ**2 

I (I , J+liK) /DELY( J)+PHI(I ,J+1 ,K)/DELY( J) )/ 

Y(J)  ) 

)*-*2-2.0E0*-U(I  ,  J,K)**-2+UU-l,  J,K)**2)/DELX*»2 
)**2/DELY( J) 
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42 


C2=(V(  I, J 
C3=V(I  ,  J+ 
C  'J=-C3 
C=2.0E0*( 
D  =  (  W(  I  »  J, 
E2A=U( 1+1 
E2B  =  U(H-1 
E2C=U<  1  +  1 
£2C=-E2C 
E2=DELY(J 
EiA=U( 1-1 
Elb=U(  1-1 
E1C  =  U<  1-1 
E1C=-E1C 
E1=DELY(J 
E  =  (E2-E1) 
F2A=V( I,  J 
F26^V(I  ,J 
F2C-V( I ,J 
F2C=-F2C 
F2=DELY(J 
F1A=-V( I, J 
F1B=V(  I, J 
F1C  =  V(  I  ,J 
F1C=-F1C 
Pl  =  OELy  (J 
F=(F2-F1) 
G2=U(  1  +  1, 
G1=U(I-1, 
G=(G2-G1) 
H=(2.0E0/ 
APHI=( 1.0 
RESI0=ABS 
EkR=MAX1( 
PHKI  ,  J,K 
CONTINUE 


,K)**2)*(1 
1 ,K)**2/DE 

C1-C2+C3)/ 
K+l)**-2-2. 
fJ+lfK)*V( 
f.JtK)  .*V(I  + 
,.J+1,KM  V( 

)*DELY(J)*' 
,J+1,K)*V( 
,.J,KH  V(  I- 
,  J+ 1  ,  K  J  *  V  i 

)*DELVCJ)* 

/(2.0E0*DE 
+lt K+i)*W< 

,  K+ 1 )  f '  Im  (  I  , 
♦  It  K+1)*W( 

)*DELY{ J)* 
* 1 ,  K  )  *  W  (  I  , 
t K)¥U( I* Jv 
+lfK)lW(I  f 

)*DELY< J)* 

/(2.0E0*DE 
J,K+1)*W( I 
Jt K+1)*W( I 

/(4.0E0*DE 
DELX**2+2  . 
EO-OMEGA)* 
(  APHI-PHK 
ERR, RES  ID) 
)  =APHI 


.0E0/0ELY(J)+1.0/DELY(  J) ) 
LY(  J) 

(DELY(  J)-t-DELY(J)) 

OEO*  W(I  f  JfK)**2  +  W(If  J,K|-**2)/DELZ**2 

I  +  lf  J+1,K)/DELY<  J)**2 

i,J,K)  Hl.OEO/DELYU  )** 2-1 .OEO/ DELY (  J)**2) 

I+ltJ+lfK)/DELY(J)**2 

(E2A  +  E2B-E2C)/(DELY(  J  )+DELY(J) ) 

I-lf  J-t-l,K)  /DELY(  J)**2 

It Jf K)*(  1.0  EC/ DEL Y( J )**2-l. OEO/DEL Y( J)** 2) 

I-lf J+1,K)/DELY(J)**2 

(E1A+E1B-E1C)/(DELY( J )+DELY(J) ) 

LX) 

It  J+1,K+1)/DELY(J1**2 

J,K+1)*(1.GE0/CELY(J  )**  2-1  .OEO/DEL Y(  J  )  +  *-2) 
I » J+1,K+1)/DELY( J)**2 

(F2A  +  F2B-F2C)/  (  DELY  (  J  )  +  DELY  ( J  )  ) 
J+1,K)/DELY( J)**2 

K)*(1.0EG/DELY(J)**2-1.0E0/DELY(  J)»<*2) 
J+l  f  K)/DELY(  jy**2 

(F1A+F1B-F1C)/(DELY( J )+DELY(J)  ) 

LZ) 

+  1,  J,K+1)+U(  I-1,J,K)*W(I-1,  J,K) 

-If  J,K+1  )+U(  I+lf.J,K)*W(  1+1,  Jf  K) 

LX*DELZ) 

0E0/DELX*^2+2.0E0/(DELY(J)*DELY(J))  ) 

PHK If J,K)+OMEGA*MA+B+C  +  C+2.0E0*(E+F  +  G) )/H 

I  ,  J , K)  ) 


C      CALCULATION  OF  PRESSURE  FROM  3D  POISSON* S  EQUATION 

C      CALCULATION  ON  A  CORNER 

C      J-l  AND  K+l  TERMS  HAVE  BEEN  ADJUSTED 

DC    43    I=2fLEl 

J  =  l 

K  =  NE 

Al=(PHi (1+1, J,K)+PHI( 1-1, J,K)) /DELX**2 

A2=(PHI (I, J,K)+PHI (I, J,K-1 ) )/DBLZ*«2 

A  3=2.  OEO*  (  PHI  (I  ,  J+1,K)  /  DEL  Y  (  J) +PHIU  t  J  +  l  ,  K  ) /DELY  (  J  )  )  / 
1 (DELY(J)+DELY(J ) ) 

A=A1+A2+A3 
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L 

C 

c 
c 


F 

F 
F 
F 
F 
F 
F 
F 
F 
F 
F 
G 
G 
G 
H 
A 
R 
E 
P 
42    C 


=  (UU«-It 
1=V(  I  t  J  + 
?=(V(  I, J 
3=VU  ,  J+ 

3=-C3 

=  2.OE0>  ( 

=  (  Ui  I  ,  J, 

2A*im  i+i 

2E=U<I+1 

2C^U(  1  +  1 

2C=-E2C 

2-£>eLY(J 

1A  =  U( I-i 

;e=u( i-i 

1C  =  U(  1-1 

1C=-E1C 

1=DELY(J 

-  {  E  2-  E  1 ) 
2A=V(  I  tJ 
2c=V(I  ,J 
2C  =  V(  I  ,J 
2C=-F2C 
2=DELY( J 
1A=V< I t J 

ie=v(i ,j 

1C  =  V(  I  ,J 
1C=-F1C 
1=DELY(J 
=.{F2-F1) 
2=U(I+lt 
1=U(I-1, 
=  IG2-G1) 
M2.0EG7 
PHI =( 1.0 
ES  ID=ABS 
RR=MAX1 ( 
HKI  »  J  ,K 
CNTINUE 


J^K)*>'2-2  .0E0*U(  I»  J,K  }*  *2+U(  1-1 , J , Kl**2)/ DELX**2 

l,K)*v2/DELY( J) 

r K)**2)*(1.0EQ/DELY(J J+l.O/DELYl J)) 

1 fKl**2/DELYCJ) 

C1-C2+C3  )/<DELY(  J)+DELYUn 

KJ**2-2.  0E0*W(I  ,  J,K)**2+W(I  ,J,K-1  )**2 )/DELZ**2 

,  J+1TK)*V(  1  +  1,  J+1,K)  /DELY(  J)^2 

».J-,K)*»V(  I+lt  JtK)»  (t.0E0/DELY(  J  ^*2-l.OE0/DELY(  J)*«2) 

iJ+l.KHVU+UJ+ltKI/DELYU)*^ 

)"DELYtJ)^ (E2A+E2S-E2C )/ (DELY( J )+CELY(J) ) 

,  J  +  1,KJ  *V(  1-1  ,J+l,K)/[)ELYf  J)*?  2 

, J,  KKV( I-li J,K)v(  I.0E0/DELYC J  )** 2-1 . 0E0/0EL Y( J ) 

t  J+lffO  +  VU-l  ,J+1,K)/DELY(J)**2 

)«DELY<J)*  IE1A+E1B-E1C)/(DELY(J)+DELY(J  )  ) 
/(2.0E0*DELX) 

+  lt  M"  W(  1  tJ+l»lO/DcLY(  J)** 2 

,K)*W(I , JtKJ*(1.0E0/DELY(J)**2-l-0E0/DELY(J )**21 

+  lf  K)*W(  I  ,J+1»K)/DFLY(  JJw-i  2 

)*>DELY(  J)*  (F2A+F2B-C2C)/(DELY( J)+CELY(J  )  ) 
*lfK-l)»N(It  J  +  lfK-l)/DELY(  JJ»*2 

fK-l  )*W(  I,  J,K^1)*(  l.QcO/DELYU I* * 2~1 . OE 0/DEL Y( J)***2) 
+ 1 , K-l ) * W ( I , J+ 1 , K-l ) / DE  L Y  U ) ♦•  *  2 


**2) 


)*DELY 
/{2.0E 
JtM*W 

JfKJvW 
/(4.0E 
DtLX** 

EO-OME 
(APHI- 


(J}< 

0*DE 
(  1+1 
(  1-1 

0*DE 
2+2. 

A) 


PHK 

ERR,&ESID) 
)  =APHI 


(F1A+F1B-F1C)/  (DELY(  J  )+DELY(J)  ) 

LZ* 

f J,K)+U(I-1» J,K-1)*W( I-1,J,K-1) 

,J,K)+U( 1+1 t JtK-l)*W ( I +1, J, K-l) 

LX*-DELZ) 

0EC/DELX**2+2.0E0/(DELY(J)*DfcLY( J)i ) 
PHI  (I  ,  J,K)+GMEGA*  (A+B+C+C+2  .0E0-'  (E  +  F  +  G)  )/H 
It J»K) ) 


CALCU 
CALCU 
J+l  A 
DO  44 
J  =  ME 
K  =  l 
A1=(P 
A2  =  (P 
A3=2. 


LA  TI ON    OF    PRESSURE    FROM    30    POISSON'S    ECUATION 

LATION    ON    A    CORNER 
NO    K-l    TERMS    HAVE    BEEN    ADJUSTED 
1=2, LEI 

HI(H-l,J,K)+PHl{I-l,J,K))/DELX*-*2 

HI (  I. J,K  +  1 )+PHI( I, JtK) )/DELZ**2 

OE  OM PHI (I  , JtK) /DELY( J)+PHI  (I ,J-1 , K)/DELY(J-1 ) )/ 
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44 


1(DELY{  J)+DELY(  J 

A=A1+A2+A3 
B=(U(I+1, J,K)** 
C2  =  (  V(  I  ,J,K)**2 
C3=V(I, J-1,K)** 
C=2.0E0*(C1-C2+ 
D=(  W(  I  ,  J, K+l)** 
E2A=U(I+1, JTK)4 
E2B-U(I+1, J,K)* 
E2C=U( 1  +  1, J-1,K 
E2=DELY(J)*DELY 
E1A  =  U(  1-1, J,K)* 
E1B=U(  1-1,  J,KH 
ElC=U(I-lf J-lfK 
E1  =  DELY(  J);J'DELY 
E=(F2-E1)/(2.0E 
F2A  =  V(  I  ,J,K+1)* 
F2B=V(  I, J, K+l)* 
F2C=V< I ,J^i,K+l 
F2  =  DELY(  JKDELY 
F1A=V{  I, J, K)*W( 
F1B=V(I  ,J,K)*-W( 
F1C  =  V(  I  ,J-1,K)* 
F1  =  DELY(J)"'DELY 
F=(F2-F1)/(2.0E 
G2=U(  1+1,  J,K+1) 
G1=U( I-i, J,K+1) 
G=(G2-G1)/(4.0E 
h=-(2.0E0/DELX^* 
APHl=(l .OEO-OME 
RESID=ABS( APHI- 
ERR=MAX1( ERRt RE 
PHI (I ,  J,K) =APHI 
CONTINUE 


-1)  ) 

Z-Z  .OEO#U( I 
l*(1.0EO/DE 
2/DELY( J-l) 
C3)/(DELY(J 
2-2.0E0*W(I 
V(  1  +  1,  JfK)/ 
V(  1+1, J,K)* 
)*V(  1+1  ,J-1 
{ J-1)*(E2A+ 
V(I-1» J,K)/ 
Mi  I-1,J,K)* 
)*V( 1-1, J-l 
I J-l)*  (E1A+ 
0*DELX) 
W( I  ,J,K+1)/ 
W( I t Jt K+l)* 
f*W(I,J-l,K 
i  J-l)^ (F2A+ 
I, J,K)/DELY 
I , J,K)*(1.0 
W(I  ,J-1,K) / 
(J-l)*(FlAt- 
0*DELZ) 
*W<  I  +  if  JfK+ 
*W(  I-lf J,K  + 
0*DELX*DELZ 
2+2.0EG/DEL 
GA)*PHI (  If  J 
PHI (I ,  J,K)  ) 
SID) 


,  J,  K  )**2+U(  I-lt  JfK)*-*  2)  /DELX*->2 
LY(J-1)+1.0/DELY(J)  ) 

)+DELY(J-l) ) 

,  J  ,  K )  *  *  2  +  W  (  I  ,  J  ,  K  )  *  *  2  )  /  D  E  L  L  *  *  2 

DELY(  J)*t2 

( 1. 0 EO/ DELY(J-l)** 2- 1.0  EG/ DEL Y(J )S*2) 

,K)/DELY(  J-l)*»2 

E2B-E2C)/(DELY( J)+DELY( J-l)) 

DELYi J)**2 

( i.OEO/DELYC J-l)*«2-1.0cO/CELY<J)**2) 

,K)/GELY(J-l)*-*2 

E1B-E1C)/(DELY( J ) +DELY ( J-l) ) 

DbLY< J)**2 

( i.OEO/DELY( J- 1)**2-1. OEO/ DELYI J )**2) 

+  l)/DELY(J-m->2 

F2B-F2C)/(DELy ( J ) +DELY ( J-l ) ) 

(  J)^2 

E0/DELY(J-1)**2-1.0E0/DELY( J)**2) 

DELY( J-l)**2 

F1B-F1C)/(DELY(  J)+DELY(  J-D) 

1)  +  U(  I-1,J,K)*-  W(I-if  J,K) 
1)+U(  I  +  lfJf  K)-*W(  I+lf  Jf  K) 
) 

X**2+2.0E0/(DELY(  J)«  DELY(  J-l)  )) 

,  K)+0MEGA*(A+B+C  +  D+2.0E0-HE  +  F  +  G)  )  /H 


CALCULATION    OF    PRESSURE    FROM    3D    POISSON'S     EQUATION 
CALCULATION    ON    A    CORNER 
J+l    AND    K+l    TERMS    HAVE    BEEN    ADJUSTED 
DO    45    I =2, LEI 
J  =  ME 
K  =  NE 

A1=(PHI ( I+ltJ,K)+PHI( I-1,J,K)) /DELX**2 
A2=(PHI  (  I,  J,K)+PHI  (  I,  J,K-1)  )/DELZ*t2 

A3=2.0E0*(PHI (I , J,K)/DELY( J)+PHI  (  I ,  J-l , K )/ DELY ( J-l )  )/ 
HCELYU  )+DELY(  J-l)  ) 
A=A1+A2+A3 

B=(  U(  I  +  lt  JfK)  **  2-2- OEO*' IM  I  f  J  f  K)  **2+U  ( 1-1  f  J,K)**2  )/DELX*=*2 
C1=V(I, J,K)**2/DELY(J) 
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45 


46 


47 


813 


C2=(V( I 

C3=V( I, 

C=2.0E0 

D=(W( I, 

£2A=U(I 

E2B=U( I 

E2C=U( I 

E2=DELY 

E1A=U( I 

E1B=U(I 

E1C  =  U(  I 

E1=DELY 

E=(E2-E 

F2A  =  V(  I 

F26=V(  I 

F2C=V(I 

F2=DELY 

F1A=V(  I 

F1B=V(I 

F1C=V(  I 

F1=DELY 

F=(F2-F 

G2=U(  1  + 

G1=U(I- 

G=(G2-G 

H=(2.0E 

APHI=(1 

RESID=A 

E^R=MAX 

PHKI  T  J 

CONTINU 


J-1,K)**2/ 
*(C1-C2+C3 
JfK)**2-2. 
+1, J,K)*V( 
+1, J,K)*V( 
+1, J-l.K)* 
(J)*DELY( J 
-1,.J,K)*V( 
-1, J,K)*V( 
- 1 ,  J- 1  »  K )  * 
(J)*DELY( J 
1)/(2.0E0* 
f  Jf.K)*W(  I  , 
iJtK)*rt( It 
, J-1,K)*W( 
(J)*DELY(  J 
, J,K-1 )~W { 
,  J  ,  K- 1 }  *  W  ( 
,J-1,K-U* 
(J)*DELY( J 
1)/(2.0E0-* 
If  J,K)*W(  I 
1,J,K)*W( I 
l)/(4.0EO* 
0/DELX**2+ 
.OEO-OMEGA 
BS( APHI-PH 
KERR,  RES  I 
,K) =APHI 


(  l.OE 
DELY( 
)  /(DE 
OEO*W 
1+1,  J 
I  +  lt  J 
V(  1+1 
-D*  ( 
I-lf  J 
1-1,  J 
V(I-1 

DELX) 
J,K)  / 
JfK)* 
I  ,J-1 
-1)M 
If  JfK 
IfJfK 
W(  It  J 

-1)*( 

DELZ) 
+  ltJt 
-ltJt 

QELX* 
2.0E0 
)*PHI 
Id, J 

D)      ^ 


0/DELY(  J-1)+1.0/DELY(  J)  J 

J-l) 

LY(  J)+DELY(J-1)) 

(If Jt KJ**2+W( U JfK-l)**2)/DELZ*'2 

,K)/DELY(  J  )*-*2 

,«)*•  (1.0EO/DELY(  J-l  )**2-1.0E0/CcLY(Jl*4  2) 

LJ-ltK)/DELY(  J-l)**  2 

E2A+E23-E2C)/(DELY(J  )+DELY(J-l)  ) 

,K)/DELY( J)*+2 

,  K)M  1.0E0/DELY(J-1)**2-1.  OEO/DEL  Y(  J  )**21 

,J-1 ,K)/DELY(J-1)**2 

E1A  +  E  1B-E 1C)/(DELY( J)+DELY( J-l)) 

DELY( J)**2 

(  1.0E0/DELYU-1)*'*2-1.0E0/DELY(J>?*2) 

,K)/OELY( J-l)**2 

F2A+F26-F2C)/(DELY( J)+DELY( J-l)) 

-1)/DELY(  J)*v-2 

-1  )*(1.0E0/DELY(J-1)**2-1.0E0/DELY(J  )-=*2) 

-1,K-1)/DELY(  J-l)**2 

F1A  +  F18-F1C)/(DELY(J)+DELY(  J-l)) 

K)+U(I-1,J,K-1)*W(I-1,J,K-1) 
K)+U( 1+1, J,K-1»*W( I+1,J,K-1) 

DELZ) 

/DELX*^2+2.0E0/(DELYU)*DELY(J-1) )) 
( It JtKi *  OMEGA* (A+B  +  C+C+2.0EO>(E  +  F  +  G)  )/H 
,K)  ) 


DO  46  J=1,ME 
DC  46  K=1,NE 
PHI  (1,  J,K)=3.0E0*PHI(2, J , K )-3 .0 EO^PHI (3 ,  J  ,  K )  +  PHI < 4, J ,K  ) 

CONTINUE 

DO  47  J=1,ME 
DO  47  K=1,NE 
PHI  (LE,  J,K)=3.0E0*-PHI  (  LEI ,  J  ,  K  )  -3  .0E0*P  H  I  ( L  E2,  J  ,  K  )  +PH  I(  L  E  3,  J ,  K  ) 

CONTINUE 

IP(ERR.GT.ERRTOL)  GO  TO  200 

IF(IC.EQ.O)  GO  TO  150 

WRITEI 6,813) 

FORMAT* 1H0,20X, 'OUT  OF  CALCULATION  POISSON') 

IF( ITER.EO.l )  GO  TO  49 
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IT  ER=1 

GC    TO    18 

49  T=T+DELT 
WFITE(6,51)    T 

51  FORMATt 1H0,20X,» TI ME    =«,E15.5) 
DO    52    J=l,50 

I=LE+1-J 

WRITE( 6,50)     U(LP1,J),U(LP2,J),U(J,1),U(I,1) 

50  FORMAT (  1H0,4E20.5) 

52  CONTINUE 
GC  TO  14 
END 

C      /*  CARD  MUST  APPPEAR  HERE 

//G0.FT06F001    DD    SYSQUT= A , SP ACE=( CYL , ( 5 , 1 ) ) 

//GC-.FT16F001    DD   UN IT=SY SDA, SPACE= ( C YL , ( 5, 1) , R LSE ) , 

//       DCB=(RECFM=VS,LRECL=3504,BLKS  IZE=3508 ) , DIS P  =  N EW , DSN  =  KORE AN  1 

//GC.FT16F001    DD    UN  I T= SYSDA , SPACE =(C YL , ( 5, 1)  , R LSE ) , 

//       CCB=(PECFM=VS,LRECL=3504,BLKS  IZE= 3508) , DI SP  =  N EW,DSN  =  K0REAN3 

//G0.FT19F001    DD    UNI T=SYSDA , SPACE= <CY L , (5 , 1 ) , RLS E ) , 

//       CCB=(RECFM=FB,LRECl=14,BLKSIZE=35  00) , DI SP  =  NE W  ,DSN  =  K0REAN4 

//G0.FT20F001    DD    UNIT  =  SYS DA,S PACE= ( CYL, ( 5, 1 ) ,  RLS E ), 

//       DC8=(RECFM=FB,LRECL=14,BLKSIZE=35  00I ,D ISP=NE W  ,DSN  =  K0REAN5 
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APPENDIX      E 
I.  FLOW  CHART  FOR  2-D   PROBLEM 


SET   I.  C. 

ITER=1 

4 

^ 

0 

TEST  ICALC 

1 

v 

ITER  =  0 

STORE   U 

4 

^ 

COMPUTE    Uf 

ITER= 1 

STORE    Ut 

r\ 

TEST  ITER 

U 

V 

1 

DELT=DELT   1 

DELT=DELT2 

f            * 

0 

TEST  ICALC 

1 

V 

1 

READ   U 

READ    Uf 

o    > 

COMPUTE   U 

COMPUTE  V 

COMPUTE   P 

TEST  ITER 

w 

1 

w 

rtiiTniiT 

U  U  1 

ru  i 

V 

6k 


APPENDIX    E 

II.    LISTING    GF    COMPUTER    PROGRAM    FOR    2-D    CASE 

THE    FOLLOWING    FORTRAN    CODED    FINITE    DIFFERENCE    COMPUTER    PROGRAM 
SOLVES    THE    FULL    TIME    DEPENDENT    NAVIER    STOKES    EQUATIONS    FOR    THE 
TwG    DIMENSIONAL    FLOW    ABOUT     A    FINITE,     INFINITELY    THIN    FLAT    PLATE 
IMPULSIVELY    STARTED    IN    ITS    OWN    PLANE.     "LEAP-FROG    TIME-WISE    INTEG- 
RATION    IS     INCLUDED    AS     AN    OPTION.    EXTENSIVE     DISC    WRITING    IS    EMPLOY- 
ED   TO    SAVE    COMPUTER    CORE    SPACE.     ADDITIONALLY,    THE     EQUATIONS 
THEMSELVES    ARE    APPLIED    TO    THE    BOUNDARIES    OF     THE    CALCULATION    REGION 
IN    LIEU    OF    EXTRAPOLATING    THE    VARIABLES    CALCULATED    WITHIN    THE 
REGION     TO    THE    BOUNDARIES. 

THE    FOLLOWING    PARAMETERS    MUST    BE    SPECIFIED 
ICALC    LE    "ME     LP1    LP2    RE    DELX    DELY     CMEGA    ERRTOL    PE    QE 

4  *  *  *  *  :&.  v  *  **  *  ****$  *$  :'..  *  *  M  *t  :.•  *  *  *!?  $$  $  *  *  $  *  a{  si  $  :«£;>;iii*  *  **:>:*-* 

ICALC    SPECIFES    THE    USE    OF    "LEAP-FROG"    INT      * 
*       LE    IS    THE     LENGTH    O*1    THE    REGION    IN    DELX    STEPS 


*  ME    IS    THE    HEIGHT    OF    THE    REGION     IN    DELY    STEPS 

•A 

*  LP1     IS    THE    START     OF    THE    PLATE 

*  LP  2    IS    THE     END    OF    THE    PLATE 

*  RE    IS    TiHE    REYNOLDS    NUMBER    BASED    ON   U    AND    L       * 

*  DELX     IS    THE    GRID    SPACING    IN    X 

*  DELY    IS    THE    GRID    SPACING    IN    Y 

*  fe 

*  CMEGA  IS  THE  RELAXATION  PARAMETER  IN  POISSONS  * 

*  ERRTOL  IS  THE  ERROR  TOLERANCE  IN  POISSONS 

*  '■* 

*  PE  IS  THE  FREESTREAM  PRESSURE  IN  P SF       * 
t  * 

*  QE  IS  THE  FREESTREAM  DYNAMIC  PRESSURE  IN  PS F 

-v  * 

*  NOTE    THAT    DELY     IS    AN    ARRAY    DELY(ME)  r' 
* 

*  ALL    VARIABLES    ARE    DIMENSIONED    F(LE,ME)  .* 

*  t. 

*  *  ?v  •«  *  *  -  **************  $   i  *  ■■■  if  ******  fc  f.  *******>£*****  ft 
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c 
c 
c 


1A**2+S  **2I) 


D  IS  EMS  I  CM    U(100,50),V(100,50),PHI(100,5G),DELY(50) 
T=O.OE0 

Ltl=LE-l 

LE2=LE-2 
LE3=LE-3 

ME1=ME-1 
ME2-ME-2 

ME3=ME-3 

REI  =  1.0E0/'«E 
PHIE=PE /( 2.0E0* QE) 

A=1.0E0/DELX 

6=1.0E0/DELY( 1) 
DELT=1.0E0/(A+0. iE0*B+2.0E0*PEI 

DELT1=DELT 
DELT2=DELT1/2.0E0 

IC  =  0 

INITIALIZE    THE    VARIABLES    WITHIN    THE    CALCULATION    REGION 

DO    2    1=1, LE 

DO    2    J=1,ME 

U(I ,J) =1.0E0 

V(  If  J)  =  O.OEO 

PHILI,  J  )=PHIE 

CONTINUE 

SET    PLATE    TO    ZERO    VELOCITY 
DC     3    I=LPlrLP2 
U(  I*1)  =  0.0E0 
CONTINUE 

CALCULATE    v    FROM    CONTINUITY    EQUATION 

DO    4    J=1,ME1 

DC    4    I =2, LEI 

C1=U(  I  +  lf J  +  l)+U( 1  +  1, J  ) 

C2=UfI-l,J+ll+U( I-1,J ) 

V(  I,  J+1)  =  V(  I,  J)-DELY(  J)  !  'C1-C2I  /(4.0E0*  DELX) 

CONTINUE 

SET    UPSTREAM    VALUES    FOR    V    VELOCITY 
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DG  5  J=1,ME 
V(1,J)=0.0E0 

5  CONTINUE 
C 

C      EXTRAPOLATE  DOWNSTREAM  VALUES  FOR  V  VELOCITY 
DC  6  J=l ,  ME 
V(LE,  J)  =  3.0EO*V(  LEI,  J)-3.0E0* V(LE2,J)+V(  LF3,J) 

6  CONTINUE 
C 

WRITE( 6,810) 
£10     FORMAT ( lhO,20X, • INTO     IC    POISSDNS') 

r 

C      CALCULATE  INITIAL  P  FIELD  FROM  PCISSONS  EQUATION 
GO  TO  200 


C 
C 


150  IC=i 

WRITE(6  ,811) 
611    FORMATt  1H0,20X,' OUT    OF     IC     PCISSONS1) 
C 

ITER=1 
14    CONTINUE 

IF( ICALC.EQ.O)    GO    TO     18 

ITER=0 
C 

REWIMD    16 
C  STORE    VALUES    OF    PRESENT    TIME    ON    TAPE 

WRITEC16)     U 
C 


18    CONTINUE 

REWIND     19 
C 
C  CALCULATE    TIME    DERIVATIVE    OF    U     FROM    MOMENTUM    EQUATION 

DO    19    I  =2,  LEI 

DO    19    J  =  2,ME1 

A-(U(I+l,J)**2-U(I-l,J>**2)/(2.0E0*DELXf 

B1=U(  I  ,  J  +  l)"  V(I  ,  J+1)/DELY(  J)^  2 

B2  =  U(  I,  J)tV-(  I,  J  )  *(1.0E0/DfcLY<  J-l)**2-l.OE0/DELY(  J)$*2) 

B3  =  U(I  ,  J-l) *V<  I  ,  J-l  )/DELY(  J-l  )*&2 

B  =  DELY(  J)*-UELY(  J-l)     (BH-B2-B3)  /{DELY(  J)+DELY(  J-l  )  ) 

C=(PHI( 1+1 ,J)-PHI(I-1, J) )/ <2.0E0*DELX) 

D2A  =  V(  1+1,  J+l)  /DELY(  J)^-2 

D2B=V(I+1, J)*(  l.OEO/DEL Y( J-l)-'  2-1. OEO/DE LY( J ) M 2 ) 

D2C  =  VU+1,  J-1)/DELY(J-1M  -2 

D2  =  DELY( J)  K)ELY< J-l) *  (D2A+D2B-D2 C) /( DELY < J)+DELY( J-l)) 

DlA=V(I-lt J+1)/DELY(J)**2 

C1B=V(I-1,  J)M1.0E0/DELY(J-1  )-^.  2-1.0E0/DELY(J  )>  *2) 

DlC=V(I-lf  J-1)/DELY(  J-l)»*2 
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903 
19 


Di  =  DELY(J)OELY(J-l)'   (  D1A+  D1P-D  1C )  /(  DEL  Y  C  J)  +DELY  ( J-l )  1 
D=(02-D1)/  12.0E0-  DcLX) 
E  1  =  U(I  , J+1)/0£LY( J) 

E2=U(  I,  J  )*tl.0E0/DELY<  J-l ) + 1. 0E G/DF LY(  J)  ) 
EJ=U(  It  J-l  >/DELV(J-l) 
E=2.0E0*(E1-E2+E3)/(DELY(J)+DELY (J-l) ) 

UTEE=-( A+B+C+f: I • ( D-E ) ) 
WRITF(i9,903 )    UTEE 
FOFMATI  1PEI4.7) 
CONT INUE 


CALCUL 

CO    22 
J-l 

A=(U( I 
B1=U(  I 

82=U(I 

B3=U(  I 

B3=-B3 

3=DELY 

C-IPHI 

D2A=V( 

D2B=V( 

D2C=V( 

D2C=-D 

D2=DEL 

D1A=V( 

D1B  =  V( 

D1C=V( 

D1C=-D 

D1=QEL 

D=(D2- 

EL=U(I 

E2=U(  I 

E3  =  U(  I 

E=2.0t 

UTEE=- 

WRITE( 

CON  TIN 

CALCUL 
DO  23 
J=ME 
A=(U(I 
B1=U(  I 
B2  =  U(  I 
B3  =  U(  I 
B=DEL Y 


ATE 
1=2, 

+  1,  J 
,  J+l 
,  J)* 

,J+1 

(  J)- 
(  1  +  1 
I+lt 
I  +  li 
I+lt 
2C 
Y(J  ) 
I -It 
I-lt 
I-lt 
1C 
Y(  J) 
Dl)/ 
,  J+l 
,  J  ) 
t  J+l 
OME 
(  A+B 
19,9 
UE 

ATE 

1=2, 

+  1»J 
,J)* 
t  J  )  # 
tJ-1 

(  J  )  * 


TIME 
LEI 


DERIVATIVE 


U     FROM    MOMENTUM    EQUATION 


>**2-U(I-l, J)**2)/(2.0E0*DELX) 
)>iV(If  J+1)/DELY(J  )'  *2 

V(  I  ,  J)  Ml  .OEC/DELY(  J)--  ■  2-1 .0  EG/ C  ELY  (J  )* v2) 
)"  V(  I,  J+l)  /DELY(  J)     -2 

OELY( J)' (3i+B2-33)/(DELY(J)+DELY(J) ) 

tJ  )--HI  (  1-1,  J)  )/(  2.0E0"  DbLX) 

J+l )/DELY(J)«  *2 

JK(  l.OEO/DELYC  J)i  -2-1  .  OEO/Dt  L  Y  (  J  )  *  *Z  ) 

J+l )/DELYCJ )* V2 

^DELY( J)*(D2A+D2B-D2C) /( DE LY( J ) +DELY ( J ) ) 
J+l )/DELY(J )**2 

J)  M 1.0E0/DELY( J) ** 2-1 .0E0/DELY1 J )**2 J 
J+1)/DELY(J )**2 

DELY( J)* (D1A+D1B-D1C)/(DELY( J)+DELY(J)) 
(2.0F0*DELX) 
)/DELY( J) 

( 1.0E0/DELY( J)  +  l. OEO/DELY(  J)  ) 
)/DELY(J  ) 

1-E2+E3) /(DELY( J)+DELY (J) ) 
+C+PEI     (D-E) ) 
03)    UTEE 


TIME    DERIVATIVE 
LEI 


OF    J     FROM    MOMENTUM    EQUATION 


)  ■•  *2-U(  1-1,  J  )**2  )/(2.0E0>DELX  ) 

V(  I  ,  J)  /DELY( J)# *2 

V(  I,  J  )  :  (  1.0E0/DELY(  J-  1  ):*■•  2-  1.0E0/DELY(  J)*  *2) 

)'  V(  I  t  J-l  )/DELY(  J-l 1**2 

DEL  Y(  J-l)     (B1+B2-B3)  /(DF.LYC  J)+DELY(J-1)  ) 
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23 


C  =  (PH 
D2A=V 
D2B  =  V 
D2C  =  V 
C2  =  DE 
D1A=V 
01B  =  V 
D1C=V 
01=  DE 
D=(D2 
E1=U( 
E2=U( 
E3  =  U( 
z  =  2.0 
UTEE= 
WRITE 
CONTI 


H  I 
(1  + 
(J  + 

(  1  + 

LY( 

(  I- 

(I- 

(I- 

LY( 

-Dl 

I  t  J 

I»J 

I,  J 

E0' 

-(A 

(19 

NUE 


+1,J)-PHI  (  1-1.  J))  /(2.0E0*  DELX) 

1, J)/DELY(J )$*2 

1,  J)M1.0E0/DELY(  J-l)  *;  2-1.0E0/DELY(J  )'*2  ) 

1,J-1)/DELY(  J-l)*?2 

J)  *DELY( J-l)  " (D2A+D2B-D2C)/(DEL Y( J )+DELY( J-l)) 

1,J)/DELY( J)**2 

It  J)  '(  1.0E0/0ELY(  J-l)     H  2-1.0E0/DELY(  J)^2) 

1 , J-l)/DELY(J-l)*/2 

J) 'DELY( J-l)' ( D1A+D18-D1C)/(DELY( J)+DELY( J-l) ) 

)/  (2.0E0--DELX) 

)/DELY( J) 

in  l.OEO/DELYC  J-l)-H.OcO/DELY(  J)  ) 

-1  J/DELY(  J-l  J 

(E1-E2+-E3)  /(DELY(  J)+OELY(  J-l)  ) 

+B+C+REIv(0-E ) ) 

,903  )    UTEE 


2  3 


CC    23 

1=1 

UTEE= 

WRITE 

CONTI 


J=1,ME 

O.OEO 
(19,903) 

NUE 


UTEE 


29 


DC    29 

I=LE 

UTEE= 

WRITE 

CONTI 


J=l  ,ME 

O.OEO 
(19,903) 

NUE 


UTEE 


IF(ITER.EQ.O)     DELT=DELT2 
IF( ITER .EQ.l)    DELT=DELT1 
IF( ITER.EQ.O )    GO    TO    402 
IF( ICALC.EQ.O)     GO    TO   402 


REWIND    16 
READ    VALUES 
READ(16)    U 


OF    U    FROM    TAPE 


402    CCNTINU' 


410 


REWIND    19 

READ    VALUES    OF    UTEE    FROM    TAPE    AND    COMPUTE    NEW    U 

DO    410    1=2, LEI 

DC    410     J=2,ME1 

READ(19,903)     UTEE 

U(  I,J)  =  U(  I, J  )KITEE-DELT 

CONTINUE 
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DC    413    1=2  ,LE1 
J=l 

RE£D(19»903)     UT E E 
U(  I  t  J)  =U(  I  ,J)+UTEE'  OELT 
413    CONTINUE 

DC    414    1=2, LEI 
J  =  ME 

FEAD(19,903)     UT h E 
U(  I, J)  =  U( I ,J)+UTEE-DELT 
4  14    CONTINUE 


1 


DC     419     J=1,ME 

1  =  1 

U(I tJ) =1.0EO 

CCNTINUE 


DC    415    I=LPltLP2 

J=l 

U(  It J)  =  O.OEO 

415    CCNTINUE 

DG    42  0    J=1,ME 
I  =LE 

U(  It  J  )  =  3.0E0vU(LEl,J)-3.0; 
420    CCNTINUE 


0^U(LE2,J)+U(LE3,J) 


DC-   30     I  =2,  LEI 
DO    30    J=1,ME1 

UEX2=( J(I+1, J)+U( 1+1, J+l) ) /2.0c 0 
UEX1=(U(  1-1,  J  )+U<  1-1,  J  +  l)  )  /2.0E0 
UEX  =  (UEX2-UEX1 ) / (2 .0 EO*  DELX  ) 
V(I,J+1)=V(I,J)-DELY(J)-UEX 
30    CONTINUE 

DO    34    I  =  1.LE 
V(  I,1)  =  0.0E0 

34  CCNTINUE 

DO    35    J=1,MP 
V(1,J)=0.0E0 

35  CONTINUE 

DC    36    J=1,ME 

I=LE 

V(  I»J  )  =  3.0E0-£  V(L£1,J)-3.0E0    V(  LE2,J)  +  V(  LE3tJ) 

3  6    CCNTINUE 
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kRITE(6,812) 
812    FGRMATC 1H0,20X, • INTO 

2  00    CCNTINUE 

ERR=0.0E0 


calculation  poisson*  ) 


CALCU 
CO  37 
00  37 
A1=(P 
A2=2. 
1( CELY 
A  =  A1+ 
B=(U( 
C1  =  V( 
C2=V( 
C3  =  V( 
C=2,0 
C2A=U 
D2B  =  J 
D2C  =  U 
C2=0E 
D1A=U 
C1B=J 
D1C=U 
C1=DE 
D='(02 

APHI  = 
RESID 
?«R=M 
PHKI 
3  7    CCNTI 


LATI 
1=2 
J  =2 
HI  (  I 
OEO- 
(  J)  + 
A2 

I  +  lf 
I  ,  J  + 
I, J) 
ItJ- 
c0'  ( 
(  1+1 
(  1+1 
(  1*1 
LY  (J 
(  1-1 
(1-1 
(1-1 
LY(J 
-01) 
E  0  /D 
(  1.0 
=  ABS 
AXK 
t  J)  = 
NUE 


ON    OF 

t'-El 

,ME1 

+1, J )+PHI 

(PHK  I,  J  + 

DCLY(J-l) 


PRESSURE    FROM    2D     POISSCN'S    EQUATION 


(I-1,J)I /Of LX     • 2 

1)/DELY(J  )+PriI(  I,J-1)/DELY(  J-l)  )/ 

) 


J  J** 

1);  * 

* ...  2  f- 
D«  * 
Cl-C 
,  J  +  l 
t  J)* 
•  J-l 
J-OE 
,  J+l 
t  J  )  ; 
i  J-l 
)  OE 
/(2. 
ELX* 
£0-3 
(APH 
ERRi 
APHI 


2-2.0 
2/DEL 
(1.0E 
2/DEL 
2  +  C3) 
)  ■  V(  I 
VU+1 
)  :  V  (  I 
LY(  J- 
):  VU 
V(I-1 
)-  V(  I 
LY(  J- 
0E0  "0 
*  2  «■  2  . 
M^GA  ) 
I-PHI 
RESID 


EO  <U 
Y(J) 
O/DE 
Y(J- 
/(DE 
+  1,  J 
,  J)  * 
+  1,  J 
1)  *  ( 
-It  J 
•  J)'; 
-It  J 
I)-  ( 
ELX) 
OEO/ 
*PHI 
(It  J 
) 


{  I,  J  )       2+U(  I-1,J )  *-*21/DELX*  -'2 

LY( J-1)+1.0E0/DELY( J) ) 

1) 

LY(  JJ+DELYC  J-l)  ) 

+1)/DELY{  J  )■•*    2 

(l.OEG/DELYU-l)  '     2-1  .OEO/DELY(J  ) 

-1)  /DELY{  J-l)r^2 

D2A+D2B-D2C)/(DEl Y( J )+DELY( J-ll) 

+  1)  /DELYi J)**2 

(  1.0cO/OELY(  J-l)  **  2-1.0E0/DELY(  J)  - 

-1  )/DELY(  J-l)^*2 

01AfD13-DlC)/(0ELY( J)+OELY( J-l)) 

(DELY(  J)  .  DELYU-1  )  ) 

(  I  ,  JJ+OMEGA*  (A+B+C  +  2.0E0-D)  /H 

)  ) 


■2) 


2) 


CALCJLATION  OF  PRESSURE  FROM  20  POISSON1 S  EQUATION 
DC  40  I  =2,  LEI 
J  =  ME 

A1=(PHI  (  I  +  l,  J  J+PHK  I-1,J)  )/DELX**2 

A2=2.0E0MPHI  (I,  J)/DELY(  JJ+PHK  It  J-l  )/ DELY  (  J-l )  )/ 
1(CELY( J )+DELY(J-l) ) 
A=A1+A2 

B  =  (U(  1  +  1,  J)  «**  2 -2.  OEO    U(I  ,  J)+U(I-1,J)**2  )/CELX*---2 
C1=V(  I,  J)*'?2/DELY(  J) 

C2=V(I f J) * -2*(1.0E0/DELY(J-1 )  +  l .OEO/DELY( J) ) 
C3  =  V(  I  *  J-l)        2/DELYU-l) 
C  =  2.0E0  MC1-C2+C3)/(DELY(  J)+DEL  Y(  J-l)  ) 
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D2A  = 
D2b= 
D2C- 
C2  =  D 
D1A  = 
015  = 
D1C= 
D1  =  D 
D=(D 
H=2. 
APHI 
RcSI 
E-R  = 
PHK 
CCNT 


U(  I+l 
UU+1 
U(  I+l 
ELY  (J 
U(I-1 
J(  1-1 
J(  1-1 
ELY(  J 
2-D  1) 
0E0/D 
=  (  1.0 
Q=ABS 
MAXK 
It  J  )  = 
INUF 


,J)*V( 

,  J )  *■  V  ( 

t  J- 1 )  >' 

J*DELY 

,  J)W( 
,J)iV( 
.J-l)* 
)- DELY 
/(2.0E 
ELX**2 
EO-GME 
(APhl- 
ERR,RE 
APHI 


I  +  lt  J 
I  +  ltJ 
VU+1 
(J-l) 
1-1  ,J 
I-lt  J 
Vll-l 
(  J-l) 
0-DEL 
+  2  .  0  E 
GA)  '■  P 
PHI  (  I 
SID) 


)/DFLY( J)  :      2 

)*  (1.0EO/DELY(J-1)><?2-1.0EO/DELY(  J)  ■■*■*  2) 

,  J-l) /DELY( J-l)**2 

(O2A+D28-D2C)/(0ELY( J) +DE LY( J- 1) ) 
)/0ELY( J)- >2 

)*  (  1.0EO/DELY1  J-l)  j=*  2-1.0E0/DELY(J)     '2  ) 
, J-l)/DiLY( J-l)> *2 

i   (D1A+D1B-D1C)/(DELY(J)+DELY (J-l)) 
X) 

0/ (DELYCJ) ; DELYCJ-1 ) ) 
HI  (I  ,  J)+CMEGAMA+B  +  C+2.0E0>  D)/H 
,  J  )  ) 


41 


CALCJLATI 
CC  41  I  =2 
J  =  l 

A1=(PHI (  I 
A2=2. OFO* 
1(0FLY( J)+ 
A=A1+A2 
B=(U(I+1, 
C1=V( I, J+ 
C2  =  V( It  J) 
C5=V(  I  ,  J+ 
C3=-C3 
0=2.0  EO-t  ( 
D2A=U(  I  +  l 
D2B=J( I+l 
D2C=U(I+1 
C2C=-D2C 
C2=DELY(J 
C1A=U(I-1 
D1B=U( 1-1 
D1C-U( 1-1 
C1C=-D1C 
D1=DELY(J 
D  =  (D2-D1) 
H=2.0E0/D 
APHI= ( 1.0 
RESIO=ABS 
ERR=^AX1( 
PHI (I, J)= 
CONTINUE 


ON  OF 
tLEl 

+  lt  J) 

(PHI  ( 
DELY( 

J)M2 
1 )  '-  2 
■**2*( 
1 )  v  *  2 

C1-C2 
t  J+l) 
tJ)*V 
,  J+l) 


PRESSURE    FRjM    20    PGISSO^J'S    EQUATION 


)--DEL 
,  J+l) 
,J  )*  V 
t  J+l  ) 

)  >OEL 
/(2.0 
ELX*~ 
EO-OM 
(APHI 
ERRf  * 
APHI 


+  PHI ( I-lt J) ) /DEL  XI*  2 

I, J+l )/CELY(J )+^HI( I.J+1)/DELY( J) )/ 

J)  ) 

-2  .  0  EO-  U  (  I »  J  )  +U  (  I  -1 ,  J  )**2  )  /  DELX**2 

/DELY( J) 

1 .3 EO/ DELY (J )+1.0E0/DELY( J )  ) 

/DELY(J) 

+C3)/(DELY(J )+DELY(J) ) 

f V(I+1, J+l) /DELY( J)*'  2 

(  I  +  l,  J)  M  1.0E0/DELY(  J)f-2-1.0E0/DELY(  J)*«2) 

-  V(I+1,  J+1)/DELY(J  )*-2 

Y(J )-<(D2A  +  D2B-D2C)/(DELY(  J )+DELY(J) ) 

4  V(  1-1,  J+1)/DELY(  J)***2 

(  I-1,J) v(  l.OEO/DELYC J)*« 2-1.0E0/DELY( J); ^2) 

■ V( 1-1, J+1)/DELY(J )*-2 

Y(J)"(D1A+D13-Q1C)/(DELY(  J)+DELY(  J)  ) 

EO*DELX) 

2+2.  OEO/ (DEL  Y(  J)  '  DELY(  J)) 

EGA)'*PHI(I,  J)+OMEGA»  (  A  +  B  +  C +2.  OE  0;  D)  /H 

-PHI ( It  J) ) 

ESID) 


DO    46    J=1,ME 
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1  =  1 

PHI (I  , J)=3.0EO  PHI (2, J) -3.0  EG -PHI (3, J) +PH 1(4, J ) 

46  CONTINUE 

DC    47    J=1,ME 
I  =  LE 

PHI(I,J)=3.0E0-PHI(LE1,J)-3.0E0-PHI(LE2,J)+PHI(LE3,J) 

47  CCNTINUE 

I-(ERR.GT.ERRTOL)     GO    TO    200 
IF( IC .EQ.O)    GO    TO    150 


WPITE<6,813) 
813    FCRMAH  lH0f20X,«  OUT    OP    CALCULATION    POISSQN') 


IP( ITER  .EQ.l ) 
I TER=1 

GO    TO    18 


Gu 


49 


49  T=T+OELT 
WRITE(6,51  )    T 

51    FORMATC  lH0t20X,»  TIME    =',E15.5) 

DO  52  J=l ,50 

I=LE+1-J 

WRITE (6, 50)  U(LP1,J),U(LP2,J),U(J,1),U(I,1) 

50  FORMAT! 1H0,4E20.5) 
5  2  CONTINUE 


GO    TO    14 

END 

/       CARD    MUST    APPPEAR    HERE 
//G3.FT06F001    DD    S YSOUT  =  A  , SPACE  =  ( CY L ,  (5  , 1 ) ) 
//GO.FT16F001    DD    UN IT= SYSDA, SPACE= ( C YL , ( 5, 1) , R L SE ) , 
//       CC6=(REC-M=VS,LRECL=3  504,BLKS  IZE= 3  508 ) , 01 S P  =  N EW , DSN  =  KOP E AN  1 
//G0.PT19F001    DD    UNI  T=  SY  ScJ  4,  SPACE  =  (  CYL  ,  (  5  » 1 )  »  R  LS  E  ) , 
//        XBMRECFM=FB,LRECL=  14,  6LKS  I  Z  E=  3500)  ,DI  SP  =NE  W  »D  SN=K0REAN4 
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